CFD-DEM Simulation of Propagation of Sound Waves in Fluid Particles Fluidised Medium

In this work, speed of sound in 2 phase mixture has been explored using CFD-DEM (Computational Fluid Dynamcis – Discrete Element Modelling). In this method volume averaged Navier Stokes, continuity and energy equations are solved for fluid. Particles are simulated as individual entities; their behaviour is captured by Newton’s laws of motion and classical contact mechanics. Particle-fluid interaction is captured using drag laws given in literature. The speed of sound in a medium depends on physical properties. It has been found experimentally that speed of sound drops significantly in 2 phase mixture of fluidised particles because of its increased density relative to gas while maintaining its compressibility. Due to the high rate of heat transfer within 2 phase medium as given in Roy et al. (1990), it has been assumed that the fluidised gas-particle medium is isothermal. The similar phenomenon has been tried to be captured using CFD-DEM numerical simulation. The disturbance is introduced and fundamental frequency in the medium is noted to measure the speed of sound for e.g. organ pipe. It has been found that speed of sound is in agreement with the relationship given in Roy et al. (1990). Their assumption that the system is isothermal also appears to be valid.


INTRODUCTION
studied sound waves by sprinkling 'Lycopodium seeds' into an oscillating column of air within a tube to identify the nodes of the standing wave.It was found that the sound waves diminished in the presence of particles and that the sound velocity changed compared to its theoretical value in air.Later, Mallock et al. (1910) studied the velocity of sound in liquid-gas mixtures such as froths.The results also showed that the speed of sound differed in a similar way to Kundt et al. (1910).Roy et al. (1990) studied the speed of sound in a gas-fluidised bed by cross-correlating pressures measured at different points in the bed and by measuring the frequency of standing wave after a disturbance had been introduced.They found that even though the average density of the gas-particle medium was greater than air, the speed of sound was significantly lower.
The velocity of sound wave in a continuous compressible medium is given by Lamb et al. (1963), (1) Where u s the speed of sound is, p is the pressure and ρ is the bulk density.For this to apply to a two phase mixture of gas and particles a number of assumptions need to be made as given by Roy et al. (1990). 1.
The particles and gas move together, 2.
The gas is compressible and obeys the ideal gas law, 3.
The particles and gas are isothermal, 4.
The particles are incompressible.
These assumptions are similar to Mallock et al. (1910), Tangren et al. (1949) and Cambell & Pitcher. (1958).Roy et al. (1990) derived an expression for the speed of sound in a homogenous two phase media as follows.If the gas behaves ideally, then Where, V g is the volume occupied by the gas, m g is its mass, p is the pressure, R is the gas constant and T is the absolute gas temperature.Differentiating Eqn.(1) with T constant gives The bulk density of the fluidised bed is given by ( Where m s is the mass of solid, m g is the mass of gas, V is the volume of gas and solid, ρ s is the density of solids, ρ g is the density of gas and is the void fraction.As changes in total volume are only due to changes of the gas, dV g = dV.Differentiating Eqn (4) gives, (5) Which can be rearranged (using Eqns.( 5) and (3)) to give (6) The speed of sound in the two phase mixture is (from Eqn (1)) then (7) The speed of sound in an ideal gas is given by where γ is the ratio of specific heat constants.Therefore the speed of sound in the two phase medium can be related to that in the pure gas by pV m RT g g = Eqn.(8) shows that the velocity of sound in solid-gas mixture can differ dramatically from that in gas alone, e.g. for air and particles at RTP, ρ s = 1000 kg m Ϫ3 , ρ g = 1.25 kg m Ϫ3 , ε = 0.4, γ = 1.4,u so = 340 ms Ϫ1 giving u s equal to 20.7 ms Ϫ1 .Roy et al. (1990) demonstrated experimentally that speed of sound in fluidised mixture of sand and air is typically 1/30 of that in air.Similar results are reported by Bi, Grace and Zhu. (1994) where the speed of sound was found to be 10 ms Ϫ1 in a fluidised mixture of air and fine particles (50 µm diameter with a density of 1580 kg/m 3 ).The approximation made that gas and particle are in isothermal state, can be justified by computing the time required for solid and gas to attain same temperature [Roy et al. 1990].The time take for the gas to equalised in the temperature with the solid is given by Eqn.(9).(9) Where h is the gas-film heat transfer coefficient, is the number of particles, c p is the specific heat capacity of the gas whose temperature is T g , solid particles temperature is T s which is taken as constant.Eqn. ( 9) can be rearranged as given in Eqn.(10).(10) From Eqn. (10), the variables shown in square bracket represents a time constant for temperature of solid particles and temperature of gas to equalise.The time constant τ s , is then given by (11) Where N u = hd/K g is the Nusselt number.By assuming that N u = 2 (lowest possible value), particle diameter d = 70 µm, specific heat capacity c p = 1045 J/KgK, density of gas ρ g = 1.28 Kg/m 3 , and the voidage (spaces between solid particle filled by gas ratio) ε = 0.4, τ s is equal to 8.4 * 10 Ϫ6 s.This indicates that the time taken for heat transfer process between solid and gas in a two phase medium is very small e.g.fluidised bed.So, it can be justified that the particles and gas are in thermal equilibrium when a disturbance (sound wave) passes by.This assumption might not be valid in larger particles because by increasing the size of particle the time constant value increases, hence increasing the time taken by the system to reach thermal equilibrium.A similar conclusion is reached by Turton et al. (1989) and Kunii and Levenspiel. (1991).In order to verify the phenomena, similar behaviour was studies by setting up an experiment as given in Roy et al. (1990).
The aim of this work is to model the speed of sound in fluidised media using the CFD-DEM simulations.The mass flux through the bed was set so that the bed was incipient fluidised i.e. (U ≈ U MF ).To generate disturbances, the fluidised bed is lifted vertically and dropped under gravity.The fluid parameters were recorded afterwards.Finally results are compared with given theoretical relationship and conclusion is drawn.

EXPERIMENTAL VERIFICATION OF SPEED OF SOUND
) ρ ε end.A simple experiment was setup to observe same behaviour as shown in Figure 1.
The experimental setup consists of a Prespex tube with internal diameter of 5 cm and external diameter of 6 cm.The 2 ends of the tube are sealed with rubber bungs.The tube is filled with Alumina Silicate particles (diameter ≈ 50 µm).These particles were fluidised by rotating the tube.As the result of rotation particles, are exposed to centrifugal force building a relative velocity between gas and particles.This cause the particles to fluidise i.e. powder gas mixture becomes free flowing and a level (horizontal) 'free surface' forms regardless of the tilt of tube.Also a significant expansion in fluidised bed was noted before and after fluidisation.Once fluidised, an impact load is induced in the bed by striking the tube on the ground.This induces the oscillations in the fluidised bed.The frequency of which was noted by making a video at 30 Hz using Sony Handycam DCR-HC14E.A .wmv clip was captured and converted into image files.The captured images were observed to measure the wave frequency.Figure 2 shows the oscillations in the bed as captured using above experimental setup.
It was found that average time period for the experiments is 0.286 s, which corresponds to a frequency of 3.50 Hz.The wavelength of the wave can be found from the height of the bed.It was noted that bed height after fluidisation was 80 cm, which corresponds to the wavelength of 240 cm as explained in Figure 3.
The speed of sound measured using above experiment is shown below: t = 0.0 s t = 0.033 s t = 0.266 s t = 0.300 s t = 0.566 s t = 0.600 s  (12) The results found via experiment were compared with the speed of sound relationship as given in Eqn.(8).The velocity of sound in solid-gas mixture was found by putting the values of variables ρ s = 3500 kg m Ϫ3 , ρ g = 1.24 kg m Ϫ3 , ε = 0.4, γ = 1.4,u so = 340 ms Ϫ1 giving u s equal to 10.51 m s Ϫ1 .The experimental and theoretical results are found in close agreement.

CFD-DEM MODELLING
The density driven formulation has been used for compressible simulations [Anderson, (1995)].In this formulation conservative form of momentum (Eqns.( 13)), continuity (Eqn.( 14)) and energy (Eqn.( 15)) equations are discretized using central difference method and the time dependent terms are discretized using the methods of lines.Ideal gas relationship is used to relate pressure and density (Eqn.( 16)).Where ρ the density of fluid is, u x ,u y ,u z are the velocities in x, y and z direction, ε is the voidage, x, y and z are space variables, i, j and k are the grid locations and t is the time, R is gas constant, T is the absolute temperature.
The 3 rd Order Adams-Bashford [Hairer et al. (1993)] explicit scheme (Eqn.(17)) has been used to advance the conservative variables, forward in time.As density variations are very small because of low velocity, double precision variables are used computationally for storing the numbers.(17) There are several factors which require attention when using a density driven formulation for CFD simulation.Firstly, if flow velocity is very small then density variations are so low that might not be captured correctly.Secondly, for stability, grid size has to be sized according to Courant Friedrichs Lewy (CFL) [Courant et al. (1967)] number (Eqn.( 18)).Where u is the velocity, ∆t is the time for each time step and ∆x is the space between each grid point.This condition states that information shouldn't travel faster than speed of sound from one cell to another.This could results in finer grid and hence making simulation computationally expensive.Thirdly, if central differencing of the convective terms is used, compressible simulations sometimes require artificial viscosity for stability and to avoid numerical divergence.
In DEM, each particle is defined as perfect sphere but the radius of each particle can be varied.The particle material is considered to be soft, which means that during a contact energy will be stored, released and dissipated, as the material deforms as shown in Figure 4. Here, this contact force is modelled by a system of springs and dashpots as shown in Eqn.(19).
Force relating displacement and velocity is given by Cundall and Strack (1979).This scheme is the same as used by Tsuji et al. (1992).The normal, F n , and tangential forces, F t , are thus, (19) where k n , γ n , δ n is the normal spring constant, damping coefficient and displacement, k t , γ t , δ t is the tangential spring constant, damping coefficient and displacement, and µ is the frictional constant, [v i, j = v j Ϫ v i ] is the relative velocity vector between the 2 particles, n is  20)) has been used as given, (20) where E * is the equivalent modulus of elasticity, R * is the equivalent radius of particles.The tangential spring constant from is given by Mindlin and Deresiewicz's (1953) as shown, (21) where G * is the equivalent modulus of rigidity, R * is the equivalent radius of particles, δ n is the normal displacement during contact.In practice this is a gross over simplification of the tangential contact force.
The normal damping coefficient (Eqn.(22)) has been introduced as given by Tsuji et al. (1992).Tangential damping coefficient (Eqn.(23)) is given by Hoef et al. ( 2004). ( where α is related to coefficient of restitution, m * is the equivalent mass and x is the displacement of particle in space ( 23) where e t is the coefficient of restitution, m * is the equivalent mass The particles area also subjected to drag forces, which were computed using the relationship developed by Beetstra et al. (2007).
(24) where V p is the volume of particle, u is the fluid velocity, v p is the particle velocity, ε is the voidage, and β Beetstra computed using relation (Eqn.(28)) There are different ways by which the fluid-solid interaction can be simulated in a 2 phase solid gas medium.One possible way is to solve the fluid dynamics equations around each particle and to specify a no slip condition at a particle fluid interface.Integration of pressure and shear forces over particle surface would then give the interaction force.This method would be accurate but very expensive when number of particles is high (above 100 approx.).A more feasible method is to capture fluid particle interaction by computing the drag on every particle using an empirical relation and then averaging it in each fluid cell (where the fluid cell size is larger than the particles size).In this later case the empirical drag force provides the required closure.Anderson and Jackson (1967) suggested a correlation for computing force acted on a particle due to fluid (Eqn.(25)).
As due to closure, there are a fixed number of particles in each cell so total force felt by the fluid due to particle is F int = nf int , that can be written in volume average form (Eqn. (26)). ( Where ∇⋅S is the term representing deviation of pressure and viscous terms and F p is the drag force exerted by the particles on fluid computed using relationship given by Beetstra et al. (2007). ( Where V cell is the volume of cell, V p is the volume of particle, N p is the number of particles in the cell, u is the fluid velocity, v p is the particle velocity and ε is the voidage and β Beetstra is given as: Where d p is the diameter of particle and µ is the kinematic viscosity.

CFD-DEM DOMAIN AND APPLIED BOUNDARY CONDITIONS
The simulated bed setup is shown below in Figure 5.The particles were allowed to move in three dimensions, but due to the narrow domain, the fluid flow was modelled as two dimensional.This was achieved by using a single fluid cell in the 3 rd dimension.Fluid was supplied uniformly through the base with a specified mass flux (achieved by specifying the rate of change of mass flux with time).No slip for the fluid was assumed on the walls.At the exit, an outflow condition was used; for the incompressible simulation, zero gradient in each component of the velocities were assumed, and for the compressible simulation a characteristic boundary condition was applied at the outlet.Simulations were performed with nominal particle diameters of 0.15 mm.
The bed is fluidised at U/UMF ≈ 1.1 by changing the rate of mass flux at the bottom boundary.Then simulation is allowed to stay at this state for time of 1.0 s in order to allow the particles to mix and create a random mixture.The disturbance is introduced by raising the bed in y-direction by 1.5 mm (10 times of particle diameter).Then bed is allowed to drop under gravity.Sequence of steps taken to introduce disturbance are illustrated in Figure 6.
In order to avoid reflection in pressure signals, characteristics boundary conditions are applied at the top boundary [Chung, (2002)].

RESULTS AND COMPARISON
Parameter set for CFD-DEM simulations are given in Table 1.The values of pressure were noted at 6 mm above the distributor after inducing the disturbance.The pressure fluctuations .
,  On the drop of the bed under gravity, the 2 phase medium undergoes natural oscillation.These oscillations correspond to natural frequency of an open ended organ pipe.The height of the bed is equivalent to quarter wavelength as illustrated in Figure 3.The speed of sound can be computed in 2 phase medium using the relationship given in Eqn.(12).The results are compared with theoretical relation given in Eqn.(8).Results are summarised in Table 2.

CONCLUSION
A good agreement has been found in CFD-DEM and theoretical results of speed of sound in 2 phase fluidised medium.The theoretical relationship has also been verified with experimental results.This also proves that assumptions made in for the derivation of the relationship are true.These results also provide the validation of given CFD-DEM model.Also damping was found in the natural behaviour of 2 phase fluidised medium.Further work is needed to explain it appropriately.
An experiment was setup to demonstrate the oscillatory motion of a fluidised bed.Roy et al. (1990) has associated oscillations in the fluidised bed with the speed of sound.Their explanation takes into account a case analogous to organ pipe with a closed end and an open

Figure 1
Figure 1 Experiment setup to study the oscillation is fluidised bed

Figure 2 Figure 3
Figure 2 The oscillations in the fluidised bed; Captured using SONY DCR14HE Handycam at 30 Hz; Fluidised bed height was found to be equal to 80 cm

Figure 4
Figure 4 Particles model of spring and dashpot for normal and tangential contact.

Figure 5
Figure 5 CFD-DEM domain of fluidised bed simulation.

Figure 6 Figure 7
Figure 6 Sequence of steps in order to generate disturbance in the fluidised bed

Table 2 :
Summary of results