Adapting Data Processing To Compare Model and Experiment Accurately: A Discrete Element Model and Magnetic Resonance Measurements of a 3D Cylindrical Fluidized Bed.

Discrete element modeling is being used increasingly to simulate flow in fluidized beds. These models require complex measurement techniques to provide validation for the approximations inherent in the model. This paper introduces the idea of modeling the experiment to ensure that the validation is accurate. Specifically, a 3D, cylindrical gas-fluidized bed was simulated using a discrete element model (DEM) for particle motion coupled with computational fluid dynamics (CFD) to describe the flow of gas. The results for time-averaged, axial velocity during bubbling fluidization were compared with those from magnetic resonance (MR) experiments made on the bed. The DEM-CFD data were postprocessed with various methods to produce time-averaged velocity maps for comparison with the MR results, including a method which closely matched the pulse sequence and data processing procedure used in the MR experiments. The DEM-CFD results processed with the MR-type time-averaging closely matched experimental MR results, validating the DEM-CFD model. Analysis of different averaging procedures confirmed that MR time-averages of dynamic systems correspond to particle-weighted averaging, rather than frame-weighted averaging, and also demonstrated that the use of Gaussian slices in MR imaging of dynamic systems is valid.

Introduction

Despite the widespread practical importance of fluidized beds, the underlying physics of these systems is still not well understood.[1] This lack of understanding is due in part to the fact that the flows in fluidized beds are difficult to measure because the beds are opaque. Consequently, nonintrusive measurement techniques such as magnetic resonance (MR) imaging, should ideally be employed. Additionally, granular matter in fluidized beds can possess gas-, liquid-, and solid-like properties,[2] making it very difficult to model accurately. To predict the behavior of fluidized beds, two characteristics must be modeled properly: (1) the contact forces between the particles, and (2) the interaction forces between the particles and the fluid. The difficulties in measuring and modeling fluidized beds potentially lead to experimental results and simulations with a combination of known and unknown inaccuracies. Thus, to cross-validate these models and measurements, it is important to compare their results in as direct a manner as possible. Here, a 3D bubbling fluidized bed has been investigated because it provides a challenging system exhibiting unsteady flow of fluid and particles with regions both dense and sparse in particle contacts.

With the drastic increase in computational power in the past few decades, computational modeling of fluidized beds has become more widespread. An accurate model is advantageous because it can provide insight into aspects such as particle mixing, particle velocity, and bubble formation and complement experimental investigations. Currently, two forms of computational models dominate the modeling of laboratory-scale fluidized beds: discrete element modeling with computational fluid dynamics (DEM-CFD) and two-fluid modeling (TFM). The main difference between DEM-CFD and TFM is that a DEM-CFD treats particles as individual objects governed by Newtonian physics, while TFM considers the particles as a continuous phase governed by continuum mechanics. The main advantage of TFM over DEM-CFD is that it can model larger beds. However, DEM-CFD can provide more detailed and accurate results since aspects of individual particles such as location and velocity can be monitored, and DEM-CFD does not require the introduction of parameters such as particle “viscosity” and “pressure”. The main imperfection in DEM-CFD is that it does not resolve fluid flow on the subparticle level, and thus requires a drag law to describe the fluid–particle interaction forces needed to close the momentum equations for both phases. Direct numerical simulation (DNS) could resolve fluid flow on a subparticle level, and thus use the no-slip boundary condition instead of a drag law to model fluid–particle interaction more accurately. However, DNS is too computationally expensive to model laboratory-scale fluidized beds. To illustrate the computational requirements of DEM-CFD, a standard desktop computer will simulate, in one day, approximately one second of the fluidization of 30 000 particles. For the same run time on the same computer, a system simulated using DNS would have to be roughly 10 times smaller in volume and number of particles, and a system simulated using TFM could be roughly 10 times bigger.[3]

A prerequisite for relying on the predictions of a discrete element model is to validate it by experiment. Several experimental techniques such as X-ray imaging, diffusive wave spectroscopy (DWS), particle image velocimetry (PIV), positron emission particle tracking (PEPT), magnetic resonance (MR), and electrical capacitance tomography (ECT) have been used to study various aspects of fluidized beds. MR is especially powerful because it can be used noninvasively to study 3D opaque systems at fine temporal and spatial resolution and does not require tracer particles. It has been used for investigation of a variety of aspects of fluidized beds[4−7] and has also proved useful in validating numerical simulations of complex flows.[8−11] However, MR measurements are produced in the Fourier domain and the effect of the MR sampling strategy, which is determined by the pulse sequence and postprocessing methods used, can be complicated when used to study dynamic systems such as fluidized beds. For example, many MR studies of fluidization have presented information on time-averaged properties, yet the method for time averaging employed by MR can be vastly different from a naive notion of time averaging. This issue will be investigated in the present work.

To date, several comparisons have been made between discrete element models and experimental measurements of phenomena in fluidized beds.[12−21] These studies have largely been conducted on 2D rectangular beds and have focused on varying DEM-CFD parameters including the drag law used,[12,13] numerical implementation of fluid–particle interaction force,[14] contact parameters,[15−19] grid size on a 2D rectangular fluid grid[12] and wall boundary conditions.[13,15,18,19] Generally, there tend to be significant discrepancies between DEM-CFD and experimental results, even after optimization of the DEM-CFD parameters. As a result, the assumption has been that the DEM-CFD is not representing, fully, the physics of the fluidized system. However, surprisingly little has been done to investigate whether the results from DEM-CFD simulations and those from experiments have been sampled and processed in the same way so as to give a direct comparison. Thus, the present investigation is concerned with comparing results from a novel 3D cylindrical discrete element model with experimental results, varying the ways in which the DEM-CFD results were sampled and processed, rather than varying the input parameters to the DEM-CFD. Accordingly, the objectives were (1) to examine how the averaging technique used in the postprocessing of DEM-CFD data affected the output, and (2) to process the results predicted by a discrete element model in a manner similar to that used in MR imaging, thereby investigating the correspondence between the averaging and slice selection techniques used in model and experiment.

Methods

Experiment

The experimental results were those of Holland et al.,[6] and the reader is referred there for a detailed explanation of the experimental arrangement. The fluidized bed was contained in an acrylic tube (i.d. 44 mm, o.d. 60 mm) and the particles were rhoeas poppy seeds (Sutton Seeds, UK). The rhoeas seeds were 1.2 mm in diameter with an apparent density of ∼900 kg m–3, giving a minimum fluidization velocity (Umf) of 0.3 m s–1 at ambient conditions. Sufficient particles were used to give a 30 mm settled bed height. The distributor was a porous glass frit 40 mm in diameter and the bed was fluidized by humidified air, with flow rate controlled by an Omega FMA5443 mass flow controller. The pressure drop across the distributor was typically between 500 and 2100 Pa, and always greater than the pressure drop across the bed (∼200 Pa).

The pulse sequence used by Holland et al.[5] was the phase-encoded velocity imaging sequence shown in Figure 1: the motion was encoded in the phase of the observed complex signal using a sine-shaped, motion-encoding gradient. A Cartesian coordinate system was employed, with the z-axis aligned along the vertical, axial, direction. The spatial resolution in the velocity measurements was 1.04 mm in the vertical (z) direction by 0.94 mm in the horizontal (x) direction (misquoted as 1.25 mm × 0.73 mm by Holland et al.[6]). Particles were selectively refocused using a soft radio frequency pulse, forming slices 5 mm thick in the y-direction for imaging in the x-z plane, shown schematically in Figure 2. These slices were achieved by manipulating the pulse sequence such that the MR signal from the particles was weighted with a Gaussian distribution based on the position of the particles in the horizontal y-direction. The Gaussian distribution was set such that the MR signal in the center of the bed in the y-direction was weighted most heavily, with the signal originating 2.5 mm from the center in the y-direction weighted 1/10th as much as that coming from the center. A diagram of the Gaussian slices used in this experiment as compared to “top hat” slices with 5 mm thickness is shown in Figure 2. The sine gradient was applied with a half period, δ, of 0.56 ms and the observation time, Δ, was 2 ms for the velocity measurements. To obtain a time-averaged image of the velocity, 48 averages of signal maps in k-space were acquired; the time between the start of successive pulse sequences, the recovery time, was 250 ms, making the total acquisition time ∼30 min.

Figure 1

Pulse sequence modeled for MR-type postprocessing of time averaged axial velocity measurements. (Reproduced with permission from ref (6). Copyright 2008 Elsevier). The flow gradients were applied in the z-direction for axial velocity measurement. The slice gradient was modeled in the y-direction to produce an image in the xz plane. The read and phase gradients were applied in the x and z directions, respectively. The phase gradient was cycled through randomly to simulate the randomness in consecutive MR imaging excitations with regard to the passage of bubbles. MR signal data was recorded during the positive portion of the read gradient with the spin echo occurring at the midpoint of the readout.

Figure 2

Plot of 5 mm y-slice in fluidized bed (left). Plot of Gaussian and top-hat signal weighting functions versus y-position in the fluidized bed (right).

Simulations

The same fluidized bed experiment was simulated computationally. The motion for each particle in the simulated fluidized bed was governed by a model developed from that of Müller et al.,[13] based on a technique originally developed by Tsuji et al.[22] This technique combines a discrete element model[23] to simulate particle motion with computational fluid dynamics of the volume-averaged Navier–Stokes equations derived by Anderson and Jackson.[24] Thus, the particles are modeled individually, with Newtonian physics and contact mechanics governing their motion, while fluid is modeled as a continuum with properties calculated at discrete locations on a fluid grid. The technique was adapted to simulate a cylindrical fluidized bed in which the computational fluid dynamics was modeled in cylindrical coordinates while the motion of the particles was modeled in rectangular coordinates. A paper presenting a complete description of the computational model has recently been submitted for publication.[25]

The normal contact force on each particle was determined using a Hertzian model, while the tangential contact force was determined using the model of Tsuji et al.,[26] in which Coulomb’s law was introduced to account for sliding. Particle and fluid motion were stepped forward explicitly in time using the third order Adams–Bashforth scheme. An explicit scheme was used to allow for pressure waves to travel through the system for future studies of pressure oscillations in fluidized beds, and the third order Adams–Bashforth scheme was used to increase simulation accuracy and stability as compared with first- and second-order time stepping techniques.

The fluid motion was modeled using computational fluid dynamics on a cylindrical grid, invoking volume-averaged Navier–Stokes equations[24] to account for portions of fluid cells being occupied by particles. To keep fluid control volumes of constant volume, a methodology was developed in which length of the control volumes in the radial direction was kept constant, yet the angle subtended by fluid control volumes decreased with the distance of the control volumes from the center. A horizontal cross section of the fluid grid used in this model is shown in Figure 3. The force of interaction between fluid and particles was modeled using the drag law developed by Beetstra et al.,[27] chosen because it matched experimental results better than other drag laws in a previous comparison between a DEM-CFD simulation and experimental results.[12]

Figure 3

Horizontal cross section of a novel CFD grid. Different colors denote different CFD cells. More CFD cells are used in the annuli further from the center, such that the fluid cells have a constant volume to ensure the same accuracy in the volume-averaged fluid equations.

As noted earlier, the fluidized bed simulated was that used experimentally by Holland et al.,[6] thereby allowing a direct comparison between model and experiment. The bubbling bed was fluidized at twice the minimum fluidization velocity (2Umf), as determined experimentally by Holland et al.[6] This was confirmed by the predictions of the DEM-CFD of the plot of pressure drop across the bed versus superficial velocity; further details of fluidization parameters are given in Table 1. The time step for the particles was chosen such that it would be 1/20th of a typical collision time between particles in this system in order to ensure each collision would be accurately simulated. The parameters for particle contacts are shown in Table 2, and were chosen on the basis of the analysis of Müller et al.[13] for modeling the poppy seeds used in MR experiments. The initial conditions and boundary conditions for the fluid are shown in Table 3. For the CFD grid, the inlet velocity in the outer annulus of fluid cells next to the distributor was set proportionally lower than the inlet velocity for the inner annuli of fluid cells to account for, as accurately as possible, the fact that, in the magnetic resonance experiments, a 40 mm diameter distributor was used in a 44 mm diameter tube, with no inlet flow occurring between radii of 20 and 22 mm. At the outlet, a 3D characteristic boundary condition[28] was used for the momentum and continuity equations, in order to ensure that pressure waves could freely exit the system without being reflected back into it. The specifications of the cylindrical fluid grid used in the model are detailed in Table 4.

Table 1

Physical Setup of the Fluidized Bed System Simulated Using DEM-CFD

parameter	value
height of settled bed	30 mm
diameter of bed	44 mm
diameter of particles	1.2 mm
density of particles	900 kg/m3
Geldart group	D
Umf	0.3 m/s
U/Umf	2.0
number of particles	30 250

Table 2

Parameters Used for Discrete Element Model of Contacts between Particles

parameter	value
coefficient of sliding friction	0.1
Young’s modulus	1.2 × 108 Pa
Poisson’s ratio	0.33
normal damping coefficient	0.02
tangential damping coefficient	0.0001
time step	1.25 × 10–6 s

Table 3

Conditions Imposed for Computational Fluid Dynamics in DEM-CFD Simulation

aspect of simulation	condition imposed
fluid-particle interaction	Beetstra et al. (2007) correlation
fluid type	compressible
time stepping scheme	explicit 3rd order Adams-Bashforth
initial pressure	105 Pa
initial velocity	0.0 m/s
inlet velocity	inner 4 annuli: 1.83 m/s
	outer annulus: 0.89 m/s
inlet voidage	0.4
boundary velocity	full-slip
temperature	298.15 K

Table 4

Sizing of Fluid Grid and Voidage Grid Used in DEM-CFD Simulation

parameter	value
fluid grid radial spacing (dr)	4.4 mm
number of inner radial control volumes	5
fluid grid vertical spacing (dz)	5 mm
number of inner vertical control volumes	23
number of angular control volumes	inner annulus: 4
	2nd annulus: 12
	3rd annulus: 20
	4th annulus: 28
	5th annulus: 36
voidage grid horizontal spacing (dx,dy)	2.2 mm
voidage grid vertical spacing (dz)	4 mm

Turbulence modeling was not included for modeling the fluid. The Reynolds number for the system was 1900 based on the bed diameter, and 80 based on the interstices between particles. At this Re, the effects of turbulence were expected to be small, due to the large inertia of the particles as well as the fact that the turbulence would be dampened by viscous forces imparted on the fluid by the particles.[29] While some work has been conducted to implement turbulence modeling in DEM-CFD models,[30−33] it is very difficult to implement turbulence models in a way which properly accounts for the effect of particles, and thus has largely been applied for modeling more dilute systems with higher Re.[30−32]

Postprocessing of the DEM-CFD Results

DEM-CFD results were analyzed in three ways to determine the time-averaged particle velocity in each voxel of a map in the x–z plane: (1) a frame-based average in which velocities were averaged, irrespective of the number of particles in the voxels, (2) a particle-weighted average, in which the velocity was averaged based on the number of particles in the voxel, and (3) a Fourier-domain average mimicking the MR-imaging acquisition process. Each of these approaches is explained below:

Frame-Based Averaging

The positions and velocities of particles were processed at instants in time separated by 10 ms intervals to determine time-averaged particle velocities using frame-based averaging. The frame-based method averaged velocities captured at different instants or “frames” in time, with equal weighting given to each frame:Here, vavg,(x, z) is the frame-based average particle velocity in the voxel centered about the point (x, z), Nframes is the number of frames used to calculate the average. The velocity vavg,(x, z) is the average particle velocity in voxel (x,z) at time k, defined bywhere n(x,z) is the number of particles in the voxel at instant k and v is the corresponding velocity of particle l in the voxel. A particle was considered to be entirely inside voxel (x,z) if it its center resided inside the limits of the voxel. The weighting function for each particle’s velocity, w, is based on the position of a particle in the y-direction, y, to account for the fact that the time-averaged velocity is only desired for a thin slice through the fluidized bed in the y-direction, rather than the entire bed. The only weighting function studied for frame-based averaging was a “top-hat” function, as shown in Figure 2. With this weighting function, particles are weighted by a function which depends on their position in the y-direction, shown in Figure 2. With the top hat function, only particles between y = −2.5 mm and y = +2.5 mm were counted; those outside the region were not, so that:where y is the y-coordinate of particle l at instant k. Furthermore, frame-based averaging is divided into two categories, “frame 1” and “frame 2″, which correspond to potential experimental procedures which read velocity signals from particles differently; these procedures could correspond to, for example, different time-averaging methods for particle imaging velocimetry (PIV) data. In the “frame 1” methodology, a frame is not counted toward the final number of frames for a particular voxel if there are no particles in the voxel at that instant in time. In “frame 2″, a frame is counted as having a velocity of zero if there are no particles in the voxel at that instant in time.

Particle-Based Averaging

In particle-based averaging, the time-averaged velocity map is weighted by the number of particles which pass through a voxel over the averaging time interval, rather than averaging over a number of time instants. Thus,With a top-hat weighting function, the difference between particle- and frame-based time-averaging can be illustrated as follows. In particle-based averaging, if a particular voxel at a particular frame in time contains two particles, it will count twice as much toward the time average as a different frame in time when the voxel only contains one particle. In frame-based averaging, the average velocity in the voxel from these two frames contributes equally toward the final average.

Two additional weighting functions were used for particle-based averaging. First, a Gaussian function was used, as shown in Figure 2. The standard deviation, σ = 1.165 mm, was chosen such that the value 2.5 mm from the center was one-tenth of that at the center:The standard deviation was varied in a later analysis to examine the sensitivity of the results to slice thickness. Second, an analysis was made for volume-based, particle-based averaging with a Gaussian weighting function:where V is the volume of particle l. Since the experimental MR images were based on raw measurements of particles weighted by both a Gaussian function and the volume of the particles, the averaging procedure of eq 6 was expected to most closely match MR averages.

Processing of DEM-CFD Results to Match the Acquisition and Processing of Raw MR Imaging Measurements

Practical complications make it impossible for any experimental technique to follow exactly any of the basic averaging procedures for obtaining a time-averaged velocity laid out in the previous section. MR imaging is expected to follow a particle-based average with Gaussian slice weighting; however, for a variety of reasons it is unclear to what extent an MR imaging measurement will differ from a pure particle-based average. Aspects such as (a) sampling Fourier coefficients at different times over a long acquisition time, (b) imperfections in gradient ramping, (c) motion encoding arising from the image encoding gradients, (d) inability to register properly a distribution of particle velocities, and (e) signal decay during individual readouts could all contribute to artifacts and deviations from a pure particle-based average. Instead, results from the DEM-CFD were processed in a manner simulating the acquisition process of an MR experiment.

The particle positions from DEM-CFD results were used to simulate MR signals from the particles throughout the pulse sequence. The simulated signal was then processed by Fourier transformation to yield a time-averaged particle velocity image, as if it were an experimental result from a phase-encoded MR velocity determination. To simulate the excitation indicated by the black bar in Figure 1, each particle was given a signal density, ρ, proportional to its volume. The signal from each particle then varied over the course of the pulse sequence depending on its position while different magnetic gradients were applied to the system, according towhere S(t) is the simulated MR signal from particle i at time t in the pulse sequence, ρ is its signal density, and ω is the frequency at which the nuclear spins inside the particle precess at time t. The frequency ω was determined by the particle’s position and the gradients in the magnetic field, according towhere γ = 2.675 × 108 rad/(s/T) is the gyromagnetic ratio of 1H nuclei in the oil in the particles, G is the three-dimensional vector describing the applied magnetic gradients and r is the position vector of the center of a particle. To reduce computational demand in postprocessing, each particle was approximated as a point located at r. Figure 1 shows the sequence of r.f. pulses and magnetic gradients used in the experimental MR imaging and modeled in the MR-type postprocessing of the DEM-CFD results. Particle positions from the DEM-CFD simulations were recorded every 10 μs for 3.38 ms to model a full MR pulse sequence. The flow-encoding gradient, Gflow, was applied to ensure that particle velocity was imaged, rather than just particle position; this gradient consisted of two sinusoidal lobes separated in time and was applied in the axial (z) direction to ensure axial velocity was imaged. The read gradient was applied in the x-direction, as an x–z plane velocity image was desired; during the 640 μs-long positive portion of the read gradient, 64 samplings of the signal from each particle were taken at different points in k space to contribute to the signal map.

The total acquisition time of the experimental measurements was 30 min. It was not possible to simulate this length of time numerically, therefore a strategy was devised to mimic the experimental acquisition as best as possible. Three approaches were used:

(1) The number of averages was reduced from 48 to 32. Experimentally, the number of averages was determined by the required signal-to-noise ratio (SNR). The simulations are essentially noise free and therefore a reduced number of averages is justified.

(2) The same measurements of particle motion were used for the positive, negative and zero flow encoding gradients.

(3) Each excitation was initiated every 10 ms instead of every 250 ms. The reduced sampling time primarily affected the phase encoding. Phase encoding was applied in the z-direction and the magnitude of Gphase determined the value of k sampled for each simulated excitation. In experimental MR imaging, the phase gradients are incremented methodically between excitations, from −32ΔGphase to +31ΔGphase over the course of 64 excitations. There is a 250 ms period between excitations in experiments, which is much longer than the average time period for a single bubble to pass through the bed (∼60 ms). Therefore, adjacent data samples of MR signal occur at random relative to the bubble motion in the fluidized bed. Since the samples occur at random, relative to bubble motion, even with methodical sampling, the final average will be representative of the average motion in the bed. In contrast, for the MR-type postprocessing of DEM-CFD predictions, each simulated pulse sequence was started every 10 ms in order to shorten total simulation time necessary. Since the simulated pulse sequences occurred at a higher frequency than bubble passage through the fluidized bed in the simulations, methodical sampling would have caused adjacent lines in k-space to be linked to the same particle velocity associated with nearby points in bubble motion, corrupting the ultimate time-average. For this reason, the phase-encoding gradients were changed randomly from −32ΔGphase to +31ΔGphase to mimic the randomness in sampling inherent in experimental MR imaging. To further mimic time averaging and improve signal-to-noise ratio (SNR), 32 full signal images in k-space were summed.

In total, four separate simulations with different initial packings of particles were used to simulate 5 s of steady bubbling each. The 20 s of steady bubbling fluidization from this ensemble of simulations was used to attain 32 images in k-space, which were ultimately converted into one time-averaged velocity image in real space.

The slice-selective pulse and gradient were modeled in the DEM-CFD postprocessing by weighting the signal of each particle by a Gaussian function (eq 5), based on the y-position of each particle at the time of the slice-selective pulse, as shown in Figure 2. The magnitudes and durations of the gradients were set to match those used experimentally by Holland et al.[6]

After the signal map, Sunweighted, from all 32 k-space signal images was obtained, it was weighted to account for T2* relaxation around the time of the spin echo. Since the spin echoes formed experimentally at the center of k-space in the x-direction, this was accounted for in the MR-type DEM-CFD postprocessing by weighting the final signal map in the k direction according towhere td = 10 μs is the dwell time between points sampled during the simulated readout. The T2* value was set to 0.6 ms to match that measured experimentally for the poppy seeds used in the Holland et al.[6] experiment. The signal map was then transformed into an intensity map in real space via a 2D Fourier transform. The complex phase, ψ(x, z), of each pixel in the intensity map, Sweighted(k, k), was then extracted and converted into an average particle velocity, according towhere vavg,MR-type(x,z) is the MR-type time-averaged velocity in voxel (x,z), Δ = 2 ms is the duration between the two flow encoding gradients, δ = 0.56 ms is the length of the flow encoding gradients and Gflow,max is the magnitude of the flow encoding gradients. The final velocity map was produced by taking a regression for each pixel among three different phase maps produced using Gflow,max −0.045, 0, and +0.045 T/m to ensure that the velocity map was not skewed by minor velocity encoding arising from the imaging gradients.

Results

The experimental results for time-averaged particle velocity were first compared with the DEM-CFD results averaged using the MR protocol, to validate the model. The DEM-CFD results were then analyzed using the other averaging techniques to examine if MR-type averaging matched a particle-based average and how far this averaging deviated from other averaging techniques.

Figure 4 gives a comparison of the time-averaged velocities of particles in the axial direction. It compares the MR-type averaging procedure of the DEM-CFD results with the experimental MR imaging results from Holland et al.[6] The two velocity maps show a similar profile of upward particle velocity in the center of the bed and downward particle velocity along the sides. This result matches the expected “gulf-stream” profile seen in small fluidized beds in which single bubbles rise predominantly through the center of the bed. Particles are rapidly carried up with the bubbles through the center and then slowly return to the bottom of the bed, moving down near the walls. The difference map (Figure 4c) shows little difference between the DEM-CFD and the MR imaging results in the bulk of the bed, with the greatest difference in magnitude occurring at the top of the expanded bed. The DEM-CFD simulation underestimates the velocity measured by MR at the central region at the top of the bed, while it overestimates the MR velocity in the outer regions at the top of the bed.

Figure 4

Comparison of average axial velocity images from (a) experimental MR imaging[6] and (b) MR-type postprocessing of the DEM simulation. Image (c) gives the difference map between images a and b. The field of view for each image is 47 mm (z) by 44 mm (x) and the resolution is 1.04 mm (z) by 0.94 mm (x).

Figure 5 shows a comparison of the time-averaged velocities from the results of the DEM-CFD simulation using different time-averaging techniques. In all cases, the magnitude of the upward velocity exceeds the magnitude of the downward velocity. The velocity maps in Figure 5 panels b–d show little difference among the various weighting procedures used for particle-based averaging. Additionally, these particle-based averages yield a very similar velocity map to that of MR-type sampling, processing, and averaging, shown in Figure 4b. Figures 5e and 5f show that frame-based averaging using the “frame 1” and “frame 2” methods heighten and dampen the magnitude of the time-averaged particle velocities, respectively, compared with particle-based averaging.

Figure 5

Comparison of average axial velocity images of DEM simulation with different averaging methods: (a) MR-type postprocessing, (b) volume-weighted, Gaussian-slice particle-based average, (c) Gaussian-slice particle-based average, (d) “top-hat” slice, particle-based average, (e) “frame 1” frame-based average, and (f) “frame 2” frame-based average. The field of view for each image is 47 mm (z) by 44 mm (x) and the resolution is 1.04 mm (z) by 0.94 mm (x).

The analysis of Figures 4 and 5 shows that the magnitude of the time-averaged profile of velocity measured by experimental MR imaging is best matched by that of the MR-type analysis and particle-based averaging of DEM-CFD results. The maximum upward, and downward, time-averaged speeds measured by MR are still greater than those yielded by MR-type processing of DEM-CFD results. Additionally, the profile of upward velocities determined by experimental MR measurements, occurring around the axis of the bed, extends to a higher vertical position than that yielded by the MR-type processing of DEM predictions, suggesting that the DEM-CFD bed has expanded slightly less than the experimental bed.

Figure 6 shows the time-averaged axial particle velocity for the different averaging procedures for a vertical slice on the axis of the bed (top), a horizontal slice 20 mm above the distributor (middle), and a horizontal slice 35 mm above the distributor (bottom). These results are consistent with those shown in the velocity map, in that MR-type postprocessing and particle-based averaging have very similar velocities which are closest to the experimental MR imaging data, while “frame 1” averaging gives significantly higher, and “frame 2” averaging gives significantly lower, velocities.

Figure 6

Comparison of average axial velocity (m/s) images with different time-averaging methods along (a) a vertical slice at the x-direction centerline, (b) a horizontal slice 20 mm above the distributor, and (c) a horizontal slice 35 mm above the distributor. The different averaging methods compared are (i) MR-type postprocessing (◇), (ii) “top-hat” slice, particle-based averaging (□), (iii) “frame 1” frame-based averaging (○) and (iv) “frame 2” frame-based averaging (△). The averaging methods are shown in comparison to the experimental MR imaging results[6] (●).

Table 5 shows the spatially averaged sum of squares of the deviation of axial velocity for each averaging technique as compared to the MR-type postprocessing and the experimental MR imaging results according to the equationswhere Δvavg,MR–type2 and Δvavg,MR2 are the square of the deviation from MR-type average and experimental MR imaging results, respectively, of a particular averaging procedure, taken over a spatial domain consisting of Npixels. The velocity vavg(n) is the time-averaged axial velocity from that averaging procedure at voxel n. Also, vavg,MRI-type(n) and vavg,MRI(n) are the time-averaged velocity from the MR-type averaging of DEM-CFD results and the actual experimental MR imaging results, respectively. The domain of Npixels consisted of all pixels up to 40 mm above the distributor. Pixels above this height were not used because the bed in the simulation did not expand much beyond this height and if these pixels had been used, the deviations would have been dominated by this region of the bed. Additionally, the experimental results could be inaccurate near the top of the field of view because particles in this region can go undetected if particles enter or leave the radio frequency (rf) coil between the initial rf pulse and signal detection. Table 5 shows that the particle-averaging procedures produce very similar time-averages to the MR-type processing, while the frame-based techniques deviate significantly. MR-type processing matches the experimental MR imaging results the best, and particle-based averaging procedures also match the MR results well, while frame-based averaging procedures deviate significantly.

Table 5

Comparison of Time-Averaged Velocity Results for Experimental MR and DEM-CFD with Different Averaging Proceduresa

averaging procedure	squared velocity difference from MR-type post-processing (m2/s2)	squared velocity difference from experimental MR5 (m2/s2)
frame 1	45.9 × 10–3	34 × 10–3
frame 2	5.10 × 10–3	17 × 10–3
particle-based, top-hat slice	0.318 × 10–3	5.8 × 10–3
particle-based, Gaussian slice	0.316 × 10–3	5.8 × 10–3
particle-volume-based, Gaussian slice	0.314 × 10–3	5.7 × 10–3
MR-type postprocessing	0	5.1 × 10–3

The squared difference in time-averaged axial velocity spatially averaged over the voxels of encompassing the fluidized bed is given by [Δvavg2 = 1/Npixels ∑Npixels(vavg(n) – vavg,MR(n))2].

Discussion

Comparison of DEM-CFD Results with Experimental MR Results

The overall profile and magnitudes of upward and downward velocity for the experimental MR measurements (Figure 4a) and the MR-type postprocessing of DEM-CFD results (Figure 4b) are in good agreement, cross-validating the DEM-CFD model and the MR measurements. The slightly higher maximum upward and downward speeds as well as the larger maximum upward velocities seen in the MR measurements, as compared to the DEM-CFD results, can probably be attributed to either (1) inaccuracies in modeling the contact mechanics or (2) inaccuracies in modeling the drag force. Inaccuracies in modeling contact mechanics could arise from the fact that the spherical particles used in the model would contact each other differently to nonspherical seeds used in the experiment; however, it is not apparent why this discrepancy would cause the differences noticed in the results. Additionally, inaccuracies in modeling the contact mechanics could stem from the values of the contact parameters used in the DEM-CFD simulation. However, this seems unlikely since Müller et al.[13] have already found the effect of varying these parameters over wide ranges to be insignificant.

The case for the difference in results being attributed to problems in modeling drag is much stronger. Using spherical particles to model nonspherical seeds affects the degree of drag, because nonspherical particles have a higher ratio of surface area to volume, which could lead to higher drag and a larger bed expansion. Additionally, the drag law of Beetstra et al.,[27] used here and developed for modeling fluid-particle interaction with spherical particles in fluidized beds, has been demonstrated to underestimate drag by Kriebitszch et al.[34] The latter found that in a direct numerical simulation (DNS) of a fluidized bed, the drag force which would be calculated on particles using the drag law of Beetstra et al.[27] in a DEM-CFD simulation was on average 33% lower than that determined by the DNS simulation using the no-slip boundary condition. This lower average drag would account for the lower time-averaged upward axial velocity in the DEM-CFD simulation as compared to the MR experiment. A separate study[35] comparing DEM-CFD to DNS simulations suggested that the drag law in DEM-CFD simulations is most inaccurate around the sides of bubbles, since particles are not evenly distributed in fluid cells and particles can move in vastly different direction and at different velocities from one another. This suggestion can account for the fact that the predicted time-averaged velocity profile is very similar at the sides and the bottom of the bed, yet less accurate in the upper, central portion of the bed, where the largest bubbles pass through the system. Additionally, Kriebitszch et al.[34] found that the bed expansion in a fluidized bed modeled using DEM with the drag law of Beetstra et al.[27] was lower than that in a fluidized bed modeled using DNS. This observation is consistent with there being less bed expansion with the DEM-CFD results (Figure 4b) compared with the MR results (Figure 4a).

Some steps could be taken to improve the modeling of drag force in this DEM-CFD system. (1) Simulated particles could be made nonspherical, matching the shape of the poppy seeds used experimentally, and a drag law similar to that of Beetstra et al.[27] developed for nonspherical particles can be used. (2) A drag law could be developed to account for particles moving in different directions and at different velocities relative to their neighbors. (3) A filter for a drag law could be created and implemented to calculate drag in fluid cells which have uneven distributions of particles because they lie on the outskirts of a bubble.

Comparison of MR-Type Postprocessing with Other Averaging Procedures

The similarity between the profiles and magnitudes of the time-averaged particle velocities yielded by MR-type postprocessing of the DEM-CFD data (Figure 5a) and particle-based averaging (Figures 5b–d) demonstrates that MR imaging effectively yields a particle-based time-average of axial velocity. This similarity is expected because the magnitude of the MR signal is weighted by nuclear spin density, which is proportional to the volume of particles containing the nuclear spins. This similarity also demonstrates that the sampling and other aspects of MR imaging mentioned in section 4.3, which could lead to artifacts and imperfections in MR measurements, do not cause MR time-averages to deviate significantly from particle-based averages. Additionally, the similarity between Figure 5 panel b and panel c shows that, for fluidized systems with a small variance in volume of particles (e.g., σ/Vavg < 0.06), weighting an average by the number of particles rather than volume of particles has little consequence on a time-averaged velocity. The similarity between panels c and d in Figure 5 shows that weighting based on a Gaussian slice versus a top-hat slice also does not greatly affect time-averaged velocity. Further investigation showed that the time-averaged velocity results were insensitive to slice thickness when both Gaussian and top-hat slices were taken with thicknesses from 2 to 10 mm, indicating that particles in this range of slice thicknesses possess the same time-averaged flow field.

The vast difference in the magnitudes of the time-averaged velocity maps of the frame-based averaging (Figures 5e,f) and the MR-type postprocessing (Figure 5a) demonstrates that MR imaging could yield different time-averaged properties of dynamic systems than other experimental procedures, such as PIV. The increased magnitudes of upward and downward velocities in Figure 5e shows that particles tend to move at greater speeds when they do not have other particles in their immediate vicinity and an experimental procedure which sampled and processed data in a similar fashion to “frame 1” averaging would capture this effect. The lower magnitudes of upward and downward velocities in Figure 5f show that an experiment which interpreted a voxel with no signal from particles at a particular instant as having zero velocity would significantly reduce the magnitude of the time-averaged velocity compared with that obtained from MR imaging. In any case, these results emphasize that care must be taken when time-averaged experimental and numerical approaches are compared.

Comparison of Computational Time and Space Required for Different Averaging Techniques

Figures 5 and 6 demonstrate that a Fourier-based MR-type time-average of velocity is equivalent to a top-hat slice particle-based average. Given this conclusion, it is important to understand the computational time and space necessary for each of these averaging techniques. While the simulation time necessary for both techniques is identical, conducting an MR-type average required approximately 500 GB of data, since particle positions needed to be output every 10 μs. In contrast, conducting a simple particle-based average would only require approximately 3 GB since particle positions and velocities only need to be output every 10 ms. Additionally, with our computer, the postprocessing of the data for the MR-type average took approximately 1 day, while the postprocessing for the particle-based average took approximately 1 min. Thus, the time and space requirements of the two averaging techniques are several orders of magnitude different, demonstrating that it is desirable to use a simpler averaging technique when the simpler technique is equivalent to that used experimentally.

Conclusions

The 3D cylindrical discrete element model developed and used in this study was validated by the similarity in axial velocity profiles and magnitudes between an MR-type analysis of the DEM-CFD results and experimental MR results. The slight differences in velocity profiles and magnitudes appear to arise from problems in modeling the drag force, both in using spherical particles to model the drag on the nonspherical seeds and also in using the drag model of Beetstra et al.,[27] which Kriebitszch et al.[34] suggest might underpredict drag forces. These drag-related issues are deserving of further work.

When comparing DEM-CFD and experimental results, it is very important to have the postprocessing procedure of DEM-CFD match that of the experiment as closely as possible. Experimental MR sampling, processing, and time-averaging of data from dynamic systems corresponds closely to a particle-based time-average and is very different from a frame-based time-average. Here, the postprocessing procedure on the raw DEM-CFD data made a much larger difference on the ultimate results than the difference seen in studies varying DEM-CFD parameters, such as drag law.[13,27]

The numerical simulations were also used to test the effect of the excitation profile used in MR imaging. Simulations were analyzed using both Gaussian and top-hat slices with thicknesses of between 2 and 10 mm. Time-averaged velocity profiles for all of these simulated slice excitations were essentially identical. These results indicate that particles in this range of thicknesses have the same time-averaged flow field, and validate the use of Gaussian excitation profiles in MR.