Modeling Dense Particle Flow in Multistage and Obstructed Flow Receivers Using High Fidelity Simulations

. Particles are a leading contender for next-generation, concentrating solar power technologies, and the design of the particle receiver is critical to minimize the levelized cost of electricity. Falling particle receivers (FPRs) are a viable receiver concept, but many new designs feature complex particle obstructions that include dense discrete phase flows. This creates additional challenges for modeling as particle-to-particle interactions (i.e., collisions) and particle drag become more complex. To improve upon existing modeling strategies, a CFD-DEM simulation capability was created by coupling two independent codes: Sierra/Fuego and LAMMPS. A suitable receiver model was then defined using a traditional continuum-based model for the air and a granular model for the particle curtain. A sensitivity study was executed using this model to determine the relevance of different granular model inputs on important quantities of interest in obstructed flow FPRs: the particle velocity and curtain opacity. The study showed that the granular model inputs had little effect on the particle velocity magnitude and curtain opacity after an obstruction.


Introduction
The design of particle receivers is critical to enabling particle-based, concentrating solar power (CSP) technology in next-generation systems.In addition to simply heating particles to the requisite temperatures to enable high-temperature power cycles, the particle receiver must be carefully designed to minimize costs and particle losses while maximizing the receiver thermal efficiency.Additionally, the receiver must be durable and show resilience to changes in environmental conditions including wind and direct normal irradiance (DNI) transients from clouds.Here, the receiver thermal efficiency is defined as the fraction in concentrated solar radiation entering the receiver that is absorbed by the particles.To achieve these sometimes-conflicting goals, modeling and simulation has been leveraged heavily in recent designs of particle receivers [1,2].Models enable rapid design evaluation which minimizes the expense and number of complimentary experimental studies.Furthermore, as particle receivers scale to commercial sizes O(10 m), experimentally evaluating new designs may be cost prohibitive.
As such, developing accurate models of particle receivers is necessary for continued development of the technology.However, models of particle receivers are typically very complex and computationally expensive due to the wide array of physics relevant to assessing either the thermal performance or the transient behavior.Several key assumptions are often made in many receiver models to reduce the computational expense.For example, in a class of particle receivers referred to as falling particle receivers (FPRs), particle-to-particle collisions are often neglected under the assumption of dilute particle flows (particle volume fractions <10%) [3].Likewise, drag on the particles is assumed to be independent of neighboring particles for the same assumptions.The relevant particle physics for a FPR is highlighted graphically in Figure 1 with key physics that are often ignored or simplified highlighted in red circles.While many of these assumptions proved to be acceptable and even necessary historically, new developments in the technology for FPRs necessitate eliminating these assumptions.For reference, a FPR is a cavity-type receiver in which particles are dropped in a curtain via gravity past a beam of concentrated sunlight.For simple FPR designs, the aforementioned assumptions in the flow were reasonable due to low particle volume fractions in the falling curtain.However, modern FPRs now include more complex features in the particle flow including perforated plates or meshes [4] (obstructed flow) or multistage troughs [5] (multistage flow).Furthermore, in some designs particles may come in contact with the cavity walls [6].In such scenarios, particles transition in and out of dense particle flow regimes violating may of the previously used assumptions and significantly affecting particle trajectories and residence time in the concentrated sunlight.
In this paper, we discuss an expanded modeling approach for particle receivers like the FPR, that enables including additional relevant physics necessary to properly evaluate these next-generation designs.To achieve this we are coupling together two independent code bases: SIERRA/Fuego [7] (simply called Fuego hereafter) and Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [8].Fuego is a low-Mach number computational fluid dynamics (CFD) capability used to model fluid flow in and around the receiver, and LAMMPS is a molecular dynamics capability that allows for modeling the complex particle interactions through the discrete element method (DEM).This paper highlights recent efforts used to couple these two code bases creating a new and mature CFD-DEM capability.We then leverage the new capability to create a model suitable for new FPRs and conduct a model sensitivity study on a section of the receiver.The study specifically targets a FPR featuring a perforated plate to identify the most relevant DEM parameters that affect the particle flow characteristics.
The remainder of this paper is organized as follows.First, a new FPR model is described highlighting key new physics that will be included to enable capturing dense discrete phase flows.Next, the new coupled CFD-DEM simulation capability is briefly described to enable execution of this new model.Then, the model sensitivity study is described and executed on the obstructed flow receiver section to demonstrate the model and discern any impact of the granular model inputs on the quantities of interest (QoIs).Finally, the conclusions of this paper are summarized.

FPR Model Development
This new receiver model uses a Eulerian/Lagrangian framework to capture critical physics within next-generation FPRs.Here, discrete particles are coupled to a continuum model of the air inside the receiver through source terms in the momentum and energy conservation equations.For brevity, only a subset of the relevant equations for the complete FPR model are included here as they would be too extensive and outside of the scope of this paper.Readers are directed to [9] for a more thorough description of the complete equation set and their corresponding assumptions (specifically for the heat and radiation transport).This paper focuses solely on the momentum coupling between the particles and the fluid as it is the most pertinent to addressing the key omissions/simplifications used for particle-to-particle interactions and particle drag.
The Reynolds-averaged Navier Stokes (RANS) equations are used to model the air in and around the receiver.The momentum transport equations are as follows: ( ) where  is the air density,   is the time-averaged air velocity vector,  is the air dynamic viscosity,  is the pressure,   is gravity vector,   are external body forces, and −  ′   ′ are the Reynolds stresses.A two-equation turbulence model (e.g.k-ε model) is typically used to close the equations.
The external body forces   are computed from the cumulative particle drag and buoyancy forces within a fluid element of the continuum.Historically, particle drag has been computed assuming flow over isolated particle spheres based on the relative velocity between the particle and the fluid.While this works well in dilute particle flows, it ignores the effect of lubrication forces between particles created by the presence of an interstitial fluid which becomes more relevant as the particle volume fraction increases.In this approach, we superimpose lubrication forces with the typical drag effects as shown below.However, while these forces will affect the trajectory of the particles, the resulting forces on the fluid domain are equal and opposite and therefore are not included in   (assuming a relatively large fluid discretization).
Particle motion for a particle of diameter   is determined by: ( ) where  , is the particle velocity vector,   is the particle density, � −   � is the magnitude of the particle/fluid velocity difference,  , are pairwise particle interaction forces between particle  and  (up to  particles) discussed below,  , is the particle mass of particle , and   is the coefficient of drag computed by: ( ) As discussed above, lubrication forces are superimposed on the particle trajectory with Eq. ( 3) to account for the presence of the interstitial fluid in dense particle flows.They are implemented as pairwise interaction forces using a dissipative force  in pair-wise interactions of particles "1" and "2" as follows: where the terms preceded by   ,  ℎ , and   represent effects of "squeeze", "shear", and "pump" with the interstitial fluid defined for rigid spheres as they interact in [10]."Twisting" effects are assumed small and thus excluded.Note that   is the particle angular velocity and  is the distance between the two interacting particles.
Particle-to-particle interactions are also superimposed on the particle motion using a Hertzian granular model where the normal forces on a particle pair   are defined by: 4 where 3 where   is the Young's modulus,  is the particle radius,  is the Poisson ratio,  is the particle overlap,  is the particle vector separating the two particles, and  , is the normal damping term.As will be shown later, values in Eq. ( 5) can be tuned depending on which particles are interacting (e.g.particle-to-particle, particle-to-obstruction).The damping term is based on a "viscoelastic" model [11] defined by: where , , where  0 is the normal damping coefficient,  is the particle velocity, and  is the particle mass.Tangential interaction forces in a particle pair   are defined using a Mindlin [12] model: where   and   are the tangential stiffness and tangential friction coefficient, respectively.The tangential damping coefficient   is simply scaled by a constant   from  0 , and   represents an increment of the elastic tangential force.This Mindlin model calculates the accumulated elastic tangential force over the contact history to ensure it does not exceed a critical value.This removes the dependence of the particle overlap on the tangential force, and if it exceeds this threshold, the tangential force is simply rescaled such that   = −   0  +  , .Additionally, an increment of the elastic tangential force   is also re-scaled as the contact unloads as: Finally, a rolling pseudo-force   is applied using a spring-dashpot-slider model [13] as follows: ( ) where   ,   , and   are the rolling stiffness, the rolling damping coefficient, and rolling friction coefficient, respectively.Much like the tangential force model, this model sets an appropriate critical value for the pseudo-force which   will not exceed.Then, a torque  =    ×   , is applied to the two particles.
As described above, this model uses a new CFD-DEM capability enabled by coupling together Fuego and LAMMPS.At each fluid timestep, Fuego solves the momentum conservation equations (Eq.( 1)) and passes local fluid properties to LAMMPS for each particle to calculate the particle's new trajectory as defined above (Eqs.( 2)-( 8)).LAMMPS updates particle positions and velocities and returns those values to Fuego where the source terms can be computed for the next timestep.This method leverages a weak coupling strategy necessitating a small fluid time step.Note that particle positions and velocities often require significantly smaller time steps to properly resolve and LAMMPS may also execute hundreds or thousands of timesteps for each fluid timestep.

Model Sensitivity Study
A model sensitivity study was executed leveraging the new CFD-DEM capability to explore the most relevant parameters affecting critical QoI in the FPR model with obstructed flow features.The QoI identified for this study includes: the particle velocity magnitude (both mean and standard deviation) after the particle obstruction (i.e. a perforated plate) at increasing distances beyond the obstruction, and the particle curtain opacity at increasing distances beyond the obstruction.Here, particle curtain opacity is defined as a geometric line-of-sight calculation of the fraction of the background physically obstructed by particles.For the particle velocity, individual particle velocities are averaged across the simulation domain (a section of a perforated plate) in 1 mm vertical windows.For the particle curtain opacity, the values are computed in 20 mm vertical windows.Both values are computed from 1 to 90 mm from the top of the obstruction for a total of 89 velocity values (mean and standard deviation) and 4 curtain opacity values.This is summarized in Table 1 For this study, the selection of the distances beyond the obstruction are defined based on expectations of the capabilities of an ongoing model validation experimental campaign.The meshed computational domain for the model sensitivity study is depicted in Figure 2. The computational domain was 60 x 60 x 180 mm.A total of 51,840 hexahedral elements were used in the spatial discretization.A perforated plate spanning the length and width of the domain was modeled in this study in a staggered hole arrangement.The holes were 10 mm in diameter and the center-to-center distance was 15 mm.The perforated plate was constructed using static particles with 1 mm diameter.A key advantage of this technique is that the continuum mesh does not need to have prior knowledge of the particle obstruction, which is critical to the scalability of the model to full particle receivers.Unfortunately, the discrete nature of the perforated plate precluded including geometric parameters in this sensitivity study.An array of moving 1 mm diameter particles is initialized in the domain above the perforated plate as shown in Figure 3.A total of ~42,000 particles are initialized in a body-centered arrangement.The location of the particles in the first layer of the body-centered arrangement is randomized but is consistent for each sample in the sensitivity study.Periodic boundary conditions are used such that particles that fall through the bottom of the domain are relocated to the top of the domain.The particles are given a small initial downward velocity and allowed to approach a steady-state solution over 1 second of simulation time.After 1 second, the QoI listed in Table 1 were not found to change significantly and are computed for the study.
Depictions of how each QoI in Table 1 is calculated are shown in Figure 3. Figure 3a highlights the 1 mm window below the perforated plate over which the mean and standard deviation of the particle velocity is computed.For the curtain opacity, the instantaneous particle positions and their pixelated spherical silhouettes at the end of the simulation are projected onto a pixelated backwall.Then, the fraction of backwall pixels obstructed by the particles is computed at the end of the simulation to arrive at the curtain opacity.A visualization of this method being applied is depicted in Figure 3b.CARBO HSP particles are often used as the particle medium in particle-based CSP, and existing numerical and experimental studies using CARBO particles [14] are leveraged to inform the constants necessary for the granular model above where available.Small differences in the granular formulation used here do exist from the previously cited work, and the authors' best judgement is used to define nominal/initial values and ranges in the sensitivity study.Note that since particles were leveraged to define the obstruction/perforated plate, granular properties were defined for both the particle-to-particle interactions and the particle-toplate interactions.Additionally, thermal effects were not captured in this study, but the effect of changing air viscosity with temperature on the particles is varied to approximate values up to 600°C.
An incremental Latin hypercube sampling (LHS) method is used to vary a total of 18 variables in the study.Each variable, defined in Table 2 is allowed to vary for each sample within the prescribed ranges with a uniform distribution.An incremental LHS study is performed to ensure that convergence is observed in the number of samples used in the study to assess the importance of each input in Table 2. 128 samples were needed to observe convergence in the results (shown below).Several snapshots from a sample in the sensitivity study are shown in Figure 4 at 0.025, 0.035, 0.07875, and 0.5 s into the simulation.Likewise, the particle velocity magnitude and curtain opacity at increasing distance from the perforate plate are plotted in Figure 5.As shown in the figure, the mean particle magnitude after the obstruction follows the theoretical kinematic equation without drag suggesting that the presence of the air is not significant at short distances.This effect has been observed previously [15].Curtain opacity after the obstruction decreases rapidly from a peak value of ~0.98 shortly after the obstruction to ~0.92 between 60 to 80 mm after the obstruction.Pearson correlation coefficients are computed from the study to quantify the relationship between the relevant QoI (Table 1) and varied model inputs (Table 2).A Pearson correlation coefficient is a measure of the strength and direction of the relationship between the QoI and the model input normalized to values between -1 and 1.A subset of the results are plotted in Figure 6 with an increasing number of samples used in the study to show convergence.As shown, each plot shows reasonable convergence with an increasing number of samples up to 128.A relatively large number of samples were required as the relationship between model inputs and QoI was found to be very weak.Ultimately, none of the model inputs showed a strong relationship (usually defined by correlation coefficients < -0.4 or > 0.4) suggesting that critical QoI are not significantly dependent on the DEM model parameters.For the particle velocity, the most relevant values included friction coefficients, and for the curtain opacity, values related to the rolling model proved to be the most relevant.It should be emphasized, that while the particle velocity magnitude and curtain opacity are weakly related to the granular model, particle trajectories are expected to be more sensitive and are the subject of future work.

Conclusions
A new CFD-DEM simulation capability was developed for particle-based concentrating solar power technologies by coupling two mature, independent codes: Sierra/Fuego and LAMMPS.
Using this new capability, a Lagrangian/Eulerian model was created to more accurately capture the particle behavior in falling particle receivers featuring multistage or obstructed flow features.Specifically, this new model includes additional physics to capture particle-to-particle interactions and more complex particle drag mechanisms in regions of dense discrete phase flow.A model sensitivity study was conducted using the new model for a sample perforated plate.Results of the study showed that the model inputs for the granular model were not significant on the relevant quantities of interest including particle velocity and particle curtain opacity after the obstruction.

Figure 1 .
Figure 1.Depiction of the relevant particle physics in next-generation FPRs.Particle physics often neglected or simplified are highlighted in red circles.

Figure 2 .
Figure 2. Meshed computational domain for the model sensitivity study (left) and a depiction of a perforated plate and particle initialization in the domain (right)

Figure 3 .
Figure 3.A 1 mm window below the obstruction over which the particle velocity is evaluated (a), and the projected particle positions perpendicular to gravity to visualize the calculation of the curtain opacity (b)

Figure 6 .
Figure 6.Pearson correlation coefficients for the mean particle velocity (left), standard deviation of the particle velocity (center), and the curtain opacity (right) after the perforated plate. / ,

Table 1 .
Quantities of interest for the model sensitivity study

Table 2 .
Model inputs varied in the sensitivity study.Initial value is listed in specified units and the min./max.values are fractional increments of the initial value.