Project Topics
Get Complete Project Material File(s) Now! »
Relative importance of the two conduction modes
The relative importance of solid-solid and solid-fluid-solid modes on the amount of heat conducted between two particles is computed by considering the presence of two identical particles having a relative temperature of 10 K. Figure 2.5 displays the solid-solid, solid-fluid-solid, and total conductive heat fluxes as a function of the distance or overlap between the two particles. It shows that in a particle bed the dominant heat transfer mechanism is the solid-fluid-solid conduction. The relative importance of the solid-solid and solid-fluid-solid modes depends on the ratio between the solids and fluid heat conductivities. When this ratio is below ~100, the solid-solid conduction have a weak control on the particle bed macroscopic heat transfer as observed in experiments (Delvosalle and Vanderschuren, 1985).
Figure 2.5: Relative contributions of the solid-solid and solids-fluid-solid heat transfers. The abscissa represents the overlap distance, in percent of the particles radius. The ordinate is a dimensionless conductive heat flux. Each curve is normalized by the value of the solid-fluid-solid heat flux when the two particles start to touch. The black thick curve is the total conductive heat flux between the particles. The blue dashed curve is the solid-fluid-solid heat flux. The orange dashed curve corresponds to the solid-solid contribution.
Coupling between the phases
Momentum coupling
The total force exerted by the motion of the fluid to the solids includes both steady and unsteady terms, and may be expressed as:
where, ⃗⃗⃗⃗⃗ is the pressure force, ⃗⃗⃗⃗⃗ is the viscous drag force, ⃗⃗⃗⃗⃗⃗⃗⃗ is the virtual mass force, and ⃗⃗⃗⃗⃗ is the Basset force. The first force describes the effect of the fluid pressure gradient. The drag expresses the steady viscous forces applied by the fluid on the particles (or vice versa) because of their relative motions. The two lasts forces are unsteady terms only important in transient dynamics that may be neglected for many applications. They both depend on the relative acceleration between the particles and the fluid. The virtual mass describes the effect of the force needed to move a volume of fluid when a particle is accelerating. The Basset term expresses the effect of the variation in the size of the viscous boundary layer (distance over which the fluid flow is affected by the presence of the particle). Neglecting the unsteady forces reduces the momentum transfer force to the drag and pressure terms and reads:
where ( ) is the momentum transfer coefficient, and ⃗⃗⃗⃗⃗( ) is the fluid velocity at the position of
the particle . The pressure gradient can be decomposed in its hydrostatic and dynamic components.
In equation (2.50), the upper line corresponds to the Wen-Yu drag coefficient and is valid for particle volume fractions below 0.2. Above 0.2, the Ergun law gives the coefficient of momentum transfer between the fluid and the particles. This last equation can be split into two terms. The first corresponds to the viscous part and is given by a Kozeny-Carman relationship describing the viscous flow at low particle Reynolds numbers, . The second is the inertial term, which depends on the relative velocity between the two phases and comes from a Burke-Plummer equation, describing the fluid kinetics at high . The Wen-Yu drag model requires to estimate the drag coefficient, , for which empirical relationships exist with. The one used here is:
Because of the Newton’s third law, the drag force exerted by the fluid on the particles must be taken into account within the interphases momentum transfer term in Eq. (2.16). The different numerical representations of the phases (Eulerian and Langrangian) impose that the drag force at the particle scale must be averaged in space to the fluid scale. The fluid-solid momentum exchange term in Eq. 2.14 may be expressed as:
with being a generic kernel indicating the contribution of a particle located at a position ⃗ to a fluid grid node located at the position ⃗⃗⃗⃗⃗. indicates the number of particles affecting the fluid at the position ⃗⃗⃗⃗⃗.
Thermal coupling
The transfer of heat between the two phases occurs by either conduction or advection of the heat at the boundary between the two mediums. The efficiency of fluid-solid heat transfer is expressed by the Nusselt number, . Form the discrete particle point of view, the heat flux between the two phases reads:
Where corresponds to the Prandtl number (See section 6 for details). Eq. (2.55) is a experimentally determined correlation of heat transfer (Ranz and Marshall, 1952; Rong and Horio, 1999). As for the momentum, computing the heat fluxes from the solid to the fluid requires to sum the heat exchanged with each individual particles through the fluid volume as:
Further details about the momentum transfer and heat transfer interpolation and averaging are given in the next section.
Numerical solvers
Overview of the SIMPLE algorithm
To solve the fluid constitutive equations, MFIX-DEM uses the SIMPLE (Semi-Implicit Method Linked Equations) algorithm (Patankar, 1980). It is an iterative method based on successive corrections of fluid velocities and pressure field. It uses a staggered grid in which fluid pressure and velocities are stored at different positions in order to avoid the convergence to checkerboard pressure fields (Fig 2.6). An overview of the algorithm steps is presented here. For a detailed presentation of the discretization of the constitutive equations and algorithm operations, see Patankar, 1980 and Syamlal, 1998. For each fluid time step, the algorithm operations are:
i: Update the fluid physical properties. The Eq. (2.11) is used to compute the new densities and, in some runs, the fluid viscosity, according to the fluid temperature field from the previous time step.
ii. The velocity and pressure gradients are computed from the results of the previous iteration or time step (for the first iteration).
x. If both residuals are below a threshold, the time step is considered as having converged and the fluid velocity, pressure, and temperature fields are used to compute the solids dynamics in the DEM part. If the residuals are above the threshold and converges (i.e. is smaller than that of the previous iteration), the algorithm restarts from step 2 with the fluid properties computed during the present iteration. When any residual diverges, the fluid time increment, , is reduced and the iteration is restarted from step 1 with the fluid properties from the previous time step.
Figure 2.6: Representation of the staggered grid use to represent the fluid phase (in 2D here). The central control volume is indicated in gray. The dotted vertical and horizontal lines represent the edges of the neighbor control volumes. The black filled circles located at the centers of the cells represent the positions of the grid where the scalar quantities (volume fraction, pressure, density, or viscosity…) are stored. The vertical arrows represent the positions where the vertical fluid velocities ( ) are stored. Similarly, the horizontal arrows indicate where x fluid velocities ( ) are stored. The black filled squares correspond to the locations of the interpolation nodes (see section 2.4.6).
DEM solver
Computing the particle motion requires one to integrate in time their accelerations given by Eqs. (2.19) and (2.20) and ensuing velocities and positions. For that, a first-order Euler time integration scheme is used (Gear, 1971). The integral in time of the particles acceleration is approximated by:
The positions of the particles are updated as:
where is the DEM time step. For the particle rotation, the Euler time integration is:
This time integration scheme is usually preferred to higher order ones (e.g. the second order Adams-Bashforth scheme; Sundaram and Collins, 1996), because with the short solid time step used, ∆ , this integration scheme is sufficient to ensure the reliability of the results with lower memory consumptions (Džiugys and Peters, 2001; Garg et al., 2010). The reliability and efficiency of the DEM approach is linked with the solid time step used. The smaller this time step is, the higher is the frequency at which the particle velocities are updated, and the more accurate the simulation is. However, the solids iterations are time consuming and too small a time step can lead to unrealistically long computation time. The solid time step used in the simulations has thus to be set to the appropriate value to ensure both the reliability and efficiency of the simulation.
To ensure the stability of the simulation, the time step has to be related to the duration of the physical processes occurring on the particles. Two phenomena must be considered here: the particle contacts and the drag acceleration. For each of these forces, a characteristic duration may 46 Chapter 2: CFD-DEM model be defined. For the contact between two particles, it is its duration, defined as (here between two particles and ):
Then the stability time step related to the particle contact is:
where is the empirical contact stability coefficient, usually set to 50. This threshold is sufficient to ensure the reliability of the simulation for a purely granular simulation, or when the fluid has a relatively low viscosity. When the fluid has a high viscosity, the stability criterion related to the drag force needs to be accounted for as well. The drag stability time step is related to the particle response time, , which corresponds to the time for a particle to reach its terminal velocity with a constant acceleration. The particle response time is inversely proportional to the fluid viscosity:
where represents the density contrast between the fluid and solids. The drag maximum time step reads:
To illustrate the effect of these values on the reliability and stability of the simulations, let’s consider an isolated particle immersed in a viscous fluid flowing in the direction steadily at a velocity ( ). The particle starts at rest and accelerates because of the viscous drag imposed by the fluid. Neglecting the effect of gravity, the particle equation of motion can be reduced to:
Figure 2.7 displays the comparison between the analytical solution of the particle motion (Eq. 2.67) and the results of a first order Euler integration with different values for the stability coefficient . It shows that when < 0.5, the particle velocity diverges. When 0.5 < < 1, the particle velocity oscillates around the expected velocity between successive time steps and converge to the fluid velocity. At = 1, after one iteration the particle reaches its terminal velocity and the drag force is canceled. In this situation, the particle velocity terminal velocity is reproduced but the time to reach it is underestimated. When > 1, the evolution of the relative velocity between the particle and the surrounding fluid may be retrieved. The accuracy of the numerical solution is proportional to the stability coefficient. A good agreement between analytical and numerical solutions is generally obtained once > 25. For application to chemically evolved mush with high melt viscosity, this time step reach values of the 10-9 s, which is a strong limitation to the application of the CFD-DEM method with such a classical time integration scheme. This limitation will be lifted in section 4.2.1.
Interpolation schemes
An accurate description of the phase coupling in an Eulerian-Lagrangian framework requires one to interpolate the fluid properties at the particle locations and, conversely, the particles mean field to the fluid grid. This allows, for instance, the smooth calculation of the free fall of a single particle across several fluid cells. For that a second order accurate in space interpolation scheme is used (Garg et al., 2012). Here, we present the key points of the interpolation algorithms. For further details, readers are redirected to Garg et al., (2012 and 2007).
Table of contents :
Chapter 1 : General introduction
1.1 Motivations
1.2 The physical processes of mush reawakening
1.3 Imaging unrest events and seismic properties of eruptible magmas
1.4 Manuscript organization and scientific questions
Chapter 2 : CFD-DEM Model
2.1 Introduction
2.2 Governing equations of the fluid phase
2.2.1 Mass conservation
2.2.2 Momentum conservation
2.2.3 Energy equation
2.2.4 State equations
2.2.5 Effect of the presence of the solid phase
2.3 Governing equations of the solids
2.3.1 Constitutive equations
2.3.2 Contact Model
2.3.2.a Collisional interactions
2.3.2.b Frictional interactions
2.3.3 Heat transfer between the solids
2.3.3.a Solid-Solid conduction
2.3.3.b Solid-Fluid-Solid conduction
2.3.3.c Relative importance of the two conduction modes
2.4 Couplings between the phases
2.4.1 Momentum coupling
2.4.2 Thermal coupling
2.5 Numerical solver
2.5.1 Overview of the SIMPLE algorithm
2.5.2 DEM solver
2.5.3 Interpolation schemes
2.5.3.a Particle side
2.5.3.b Fluid side
2.5.4 Boundary conditions
2.5.5 Dimensionless numbers
Chapter 3 : Effects of lubrication on mush dynamics
3.1 Introduction
3.2 Method
3.2.1 Formulation of the BBO equation
3.2.2 CFD-DEM Model
3.3 Results
3.3.1 Grains scale
3.2.1.a Scaling of the relative importance of the forces exerted on a particle
3.2.1.b Dimensionless formulation
3.3.2 Macroscopic scale
3.2.2.a Experiment 1 : Rayleigh-Taylor instabilities
3.2.2.b Experiment 2 : Injection of a fresh magma into a mush
3.3.3 Interpretation
3.4 Discussion
3.4.1 Influence of the crystals size and shape
3.4.2 Comparison with other studies
3.4.3 Implication on magma rheology and magmatic system dynamics
3.4 Conclusion
Supplementary section S3.1
Supplementary section S3.2
Supplementary section S3.3
Supplementary section S3.4
Supplementary section S3.5
Chapter 4 : CFD-DEM modeling of recharge events within mush
4.1 Introduction
4.2 Method
4.2.1 Numerical method
4.2.2 Numerical setup and experiments
4.2.3 Dimensionless parameters
4.3 Results
4.3.1 Effect of buoyancy and viscosity
4.3.2 Injection velocity
4.3.3 Results summary
4.4 Discussion
4.4.1 Model limitations
4.4.2 Implication on mush dynamics and on the modeling of crystal-bearing magmas
4.5 Conclusion
Supplementary section S4.1
Supplementary section S4.2
Supplementary section S4.3
Chapter 5 : Numerical simulations of the mixing caused by a magma intruding a resident mush
5.1 Introduction
5.2 Method
5.1 Results
5.2 Discussion
5.1 Conclusion
Chapter 6 : The seismic properties of eruptible magmas
6.1 Introduction
6.2 Theoretical assumptions
6.3 Method
6.3.1 Conservative equations of the phases
6.3.2 Fluid-solid and solid-solid couplings
6.3.3 Calculation of the acoustical properties of the suspension
6.3.4 Magmas under consideration
6.4 Results
6.4.1 Magmas physical properties
6.4.2 Fluid-solid momentum coupling
6.4.3 Solids-solids momentum coupling
6.4.4 Wave velocities and attenuation coefficients in suspensions
6.4.4.a Compressional waves
6.4.4.b Shear waves
6.4.5 Compressional waves in suspensions
6.4.6 Shear waves in suspensions
6.4.7 Application to magmas
6.5 Discussion
6.5.1 Limitations
6.5.2 Model validation
6.5.3 Implications for magmas
6.5 Conclusion
Supplementary section S6.1
Chapter 7 : General conclusions
