Modeling of Magnetic-Field-Assisted Fluidization: Model Development and CFD Simulation of Magnetically Stabilized Fluidized Beds

Magnetic-field-assisted fluidization is starting to be considered as a viable alternative to standard fluidized beds for those operations (such as particle separations, filtration, adsorption) in which the solid phase can be made of magnetic particles or, alternatively, the fluidizing agent is a ferro-fluid; thus the fluid bed responds to the action of magnetic fields, and stabilized fluidization regimes can be generated. One of the major difficulties to be tackled is the development of a predictive model capable of estimating the stabilized-to-bubbling transition velocity for a given magnetic field or, on the other hand, the magnetic field intensity required to stabilize the bed to a quiescent condition. The fluid dynamics prediction of a stabilized bed is also a challenging task at the moment. On this basis, a very simple model for the description of MSFB was derived in this contribution starting from basic fluid dynamics and magnetodynamics equations. The model was implemented in a commercial CFD code in order to simulate the effect of the magnetic field onset on a freely bubbling fluidized bed.


Introduction
Magnetic-field-assisted fluidization (MFAF) is a technology that, in the last decades, has started to be considered as a possible alternative to standard fluidized beds for those operations in which magnetic particles or ferrofluids can be employed. The Magnetic Stabilized Fluidized Bed Reactor (MSFBR, or simply MSFB) is a fluidized bed in which the solid phase is composed of magnetic particles or, alternatively, the fluidizing agent is a ferro-fluid; the whole system is surrounded by an array of coils that generates a magnetic, usually axial, field.
A large number of applications for MSFBs have been proposed in the past. The main fields in which the MSFB has been applied includes particle separations, filtration, fluidized bed reactors, fluidized bed adsorption, biochemical applications, fluid bed chromatography. In fact, in its simplest form, the magnetic field action eliminates gas bubbling from a fluidized bed (i.e. magnetically 'stabilizes'), giving rise to smooth fluidization at higher gas flow rates than those achievable in the non-stabilized bed. An MSFB has the primary advantage of combining the low pressure drop of a fluidized bed with the bubble-free operation of a fixed bed (Rosensweig, 1978;Saxena and Shrivastava, 1991;Webb et al., 1996).
A systematization of the operating modes of MFAF has been done by Hristov (2002). The fluid flow and the magnetic field can be applied independently, hence two principal magnetization modes are possible: magnetization FIRST and magnetization LAST.
The term magnetization FIRST was introduced by Siegell (1987). This mode implies the application of the magnetic field on the fixed bed with an isotropic structure and static particle arrangement and the fluidization afterwards (Kirko and Filippov, 1960;Rosensweig, 1979;Hristov, 2002). In this case inter-particle forces play an important role. As a result of particle magnetization, the cohesive force of magnetic nature emerges. After the magnetic field application, an increase in flow rate gives rise to a deformation of the initial static bed with simultaneous orientation of particles along the field, leading to the socalled stabilized bed as suggested by Rosensweig (1979).
The bed appears like a fixed bed with anisotropic particle arrangement induced by field lines. This state has been considered as a transitional state between fixed bed and fluidization or as homogeneous fluidization without bubbles. Further increase of the flow rate leads to fluidization of the particle aggregates and then the formation of bubbles, with the hydrodynamic forces progressively overcoming the magnetic forces until eventually particle aggregates are destroyed and fully bubbling fluidization occurs.
In magnetization LAST mode, the magnetic field is applied to a fully fluidized bed. The magnetic force causes the particles to reduce their motion and form aggregates which are aligned along field lines. At high field intensity, the bed will ultimately collapse resulting in a fixed anisotropic structure of particle aggregates. As in the previous case, the regime of homogeneous fluidization of strings is observed.
Experimental investigation in magnetization LAST is scarce compared to magnetization FIRST. Hristov and Ivanova (1999) performed different experimental tests, obtaining similar regimes to those obtained in magnetization FIRST mode. In magnetization LAST mode, a similar sequence of fluidization regimes is observed, increasing the magnetic field intensity for a given value of the fluid velocity.
However, even if the technique of MSFBs is almost 40 years old, it has not been applied yet on a large scale. One of the main problems is the choice of a suitable magnetic system. Another one is the ability of predicting with acceptable accuracy the behavior of the full-scale equipment.
The recent development of mathematical modeling of particulate solids behavior, together with the increased computing power, enables researchers to simulate the behavior of fluidized powders (with and without applied magnetic fields) and to link fundamental particle properties directly to the powder behavior. In this regard, computational fluid dynamics (CFD) modeling provides a fundamental tool to support engineering design and research in multiphase systems. In general, many authors recognize that computational modeling in multiphase systems has the potential to increase process efficiency and reduce the number of scale-up steps in the design of reliable commercial plants.
Among the many aspects of research in this field, there is a substantial interest in trying to develop a model capable of estimating the regime transitions occurring in the MSFB, as well as to correctly simulate the fluid dynamics of the bed itself.

Literature review
Mathematical modeling of fluidized beds can be regarded as a multiphase flow problem. In a fluid bed, solid particles are suspended in a fluid, the former having a discrete nature while the latter is considered a continuum. Thus, fluidization is described as a dispersed particulate two-phase flow. The great variety of different flow regimes in fluidization (Kunii and Levenspiel, 1991;Geldart, 1973) outlines the complexity of multiphase flow. This in turn implies a major difficulty for the general formulation of governing equations. Complete descriptions of Eulerian-Eulerian models derivation and additional information can be found in Ishii (1975), Anderson and Jackson (1967) and Jackson (2001). The Eulerian-Eulerian modeling approach can be extended to the mathematical modeling of magnetofluidized beds. This implies the inclusion of additional magnetic body forces in the momentum balance equations, which have to be coupled with Maxwell's magneto-static equations. This can be done within the framework of existing models. However, two very important aspects have to be noted: the Eulerian-Eulerian approach considers the solid phase as a continuum fluid, hence magnetic interparticle force and magnetic torque cannot be directly considered, hence these models cannot predict typical phenomena of MSFBs such as the formation of strings and aggregates at high field intensities. However, the relevant macroscopic effects can in principle be accounted for, ultimately allowing the simulation of large-scale units.
Conversely, in order to simulate these phenomena, Eulerian-Lagrangian approaches would be required, able to directly account for the strong inter-particle forces, intense local polarization and complex structures deriving from intense magnetic field action (Rosensweig and Ciprios, 1991).
The first rigorous mathematical description of fluidization in the presence of a magnetic field was proposed by Rosensweig (1985). The model developed by Rosensweig is based on that proposed by Anderson and Jackson (1967). The latter follows a semi-fundamental approach based on the definitions of space average quantities on a volume which respects the scale separation condition. Rosensweig extends the analysis of Anderson and Jackson by also applying the averaging procedure to the magnetostatic equations. Another model for the description of magneto-fluidized beds was developed by Brandani and Astarita (1996). The model is based on Foscolo and Gibilaro's Particle Bed Model (Foscolo and Gibilaro, 1987), which has proved to give qualitative and quantitative agreement with experimental results. In this new model, magnetic forces acting on particles are included in the one-dimensional equations of change of Foscolo and Gibilaro. Within the field of magneto-stabilized fluidized beds (MSFB), very few contributions in the pertinent literature dealt with the CFD simulation of this class of equipment: Li et al. (2010) simulated the flow behavior of gas and solids in a two-dimensional magnetically assisted bubbling flu-idized bed with magnetic balls under a vertical-gradient magnetic field. In this case, the motion of the gas was simulated by computational fluid dynamics (CFD), while the particles were simulated using the discrete element method (DEM), thus allowing for the simulation of particle chain formation; thanks to DEM modeling it is also possible to assess the decrease of particle diffusion coefficient with the increase of magnetic field intensity.
The classic multi-fluid modeling approach was also used in the relevant literature: even if some phenomena can only be partially simulated (at most, their macroscopic effect), the relative computational easiness allows the simulation of large-scale equipment. Wang et al. (2002) carried out a numerical simulation (based on Two Fluid Model and Kinetic Theory of Granular Flow) and relevant experimental validation of a cylindrical magnetofluidized bed. The numerical results indicated that below nil or a weak magnetic field, the typical bubbling fluidization was obtained, while under a moderate magnetic field, the particles display stable fluidization by restraining bubble formation. The numerical results were found to be in agreement with the experimental data, reflecting the ability of CFD codes to tackle the simulation problem. Xu and Guan (2003) simulated an Air-Dense Medium Fluidized Bed in which a magnetic medium was used in order to further stabilize the fluidized bed and increase the equipment efficiency.
In the present work, a first attempt to develop and implement a simple mathematical model in a commercial CFD code has been made in order to perform 2-dimensional simulations of a magnetically stabilized fluidized bed. The commercial CFD code CFX 4.4 was adopted, which allows substantial flexibility in including additional terms in the momentum balance equation.

Mathematical model
In this section, the mathematical model to be used for the CFD simulation of magneto-fluidized beds is presented. This model is an extension of the fluidization model proposed by Brandani and Zhang (2006).
In general, the basic equations describing the momentum balance in a fluidized bed are usually in a symmetric form and always include a term for drag force, for gravity force and for the pressure gradient. However, this formulation is unable to predict the existence (experimentally observed) of the regime of homogeneous fluidization, so the model has to be suitably extended (by means of additional terms) in order to be able to predict the transition between homogeneous and bubbling fluidization.
In general, three types of additional forces can be added: fluid-particle interactions (other than the normal drag force), particle-particle interactions, and forces re-sulting from the averaging procedure to obtain continuum formulation. In any case, these additional terms lead to an elasticity of the bed or particle pressure which allows the existence of stable, smooth fluidization.
In the Particle Bed Model by Foscolo and Gibilaro (1987), subsequently revised in Gibilaro (2001), a fluidparticle force term is added to the solid momentum balance. However, in accordance to Newton's third law, a force arising from fluid-particle interactions should cancel out over summation and not be present in the overall momentum balance. In the model by Brandani and Zhang (2006), the discretized momentum balance equations are formulated in such a way that if a term is added to the particle-phase equation, a similar term appears also in the fluid-phase equation.
The basic relationships proposed by the authors are assumed as the basis for the present work. In this contribution, the momentum balance equations are extended here with the addition of a magnetic force term: Notably, the momentum balance equations include the usual drag force term (F D ), and two additional elastic terms (dependent on the volume fraction gradient) appear for each phase: the fluid dynamic elasticity (E p , E f ) as derived by Brandani and Zhang (2006), and the magnetic force terms (E mp , E mf ). The factor δ that appears in the last two equations is a characteristic length of the order of the particle's diameter which arises from the formulation of a discretized balance. These equations are obtained without any simplifying assumption except for neglecting higher-order terms when shifting from the discrete form to the differential form (Brandani and Zhang, 2006).
A more complex problem is that of calculation of the additional body forces due to the magnetic field. The basic relations linking the magnetic field, H, the induction field, B 0 , and the magnetization, M, are: Also, Maxwell's equations hold under the assumption that anywhere within the flow, a local region exists that is large compared to the size of any particle but small compared to the distance over which the field would vary significantly: Following the procedure employed by Rosensweig and Ciprios (1991), and by Brandani and Astarita (1996), the following general constitutive equations for the magnetic force can be derived: These equations are perfectly symmetrical, but are valid only for the two limiting cases of χ f = 0 (the solid is magnetizable, Eqns. 10-11) or χ p = 0 (the fluid is magnetizable, Eqns. 12-13). The corresponding expression for the magnetic forces is obtained by combining these equations with a constitutive relation for the bed susceptibility. It can be shown that, using the Clausius-Mossotti relationship (Banhegyi, 1986), the equations by Rosensweig and Ciprios (1991) are obtained:

Stability criterion
The kinematic and dynamic wave velocities can be derived following the linearization procedure of Wallis (1969) and Gibilaro (2001). Wallis' theory states that voidage perturbation velocities in the bed are bounded by the kinematic wave velocity and the dynamic wave velocity. The kinematic wave velocity is obtained when the steady state flow depends on the voidage, and in a fluidized bed it is only a function of the drag equation (Wallis, 1969). On the other hand, the dynamic wave velocity is the speed at which voidage perturbations run through an elastic or compressible medium, resembling, e.g. sonic waves in a gas. This fundamental result was applied by Brandani and Zhang (2006) to the limiting case without magnetic forces and results in: and the dynamic wave velocity: where, under the quasi-equilibrium approximation: The superficial velocity U is calculated from the bed equilibrium relationship and the drag force is expressed as Gibilaro (2001), with the Dallavalle expression for C D : The criterion for minimum bubbling can be expressed in terms of a kinematic and a dynamic wave velocity as Wallis (1969): the bed is stable when the kinematic velocity is smaller than the dynamic wave velocity. In the presence of additional magnetic forces acting on the particles, the kinematic velocity remains unchanged, while only the G term is modified to the following expression that has to be adopted in Eqn.17 For the investigated system, the prediction results are reported in Fig. 1 in terms of wave velocities (kinematic and dynamic) as a function of the fluid-phase volume fraction and magnetic field intensity.
When the magnetic field intensity is set to B 0 = 0.088 T or below, the minimum bubbling voidage ε mb is below the voidage at minimum fluidization conditions ε mf = 0.4 (typical of random particle arrangement): this is a characteristic of all Geldart Type B particles (always bubbling when fluidized).
Above B 0 = 0.088 T ε mf > ε mb , the particles behave as Geldart Type A powders, e.g. conditions exist in which homogeneous fluidization is feasible before bubbling occurs.
For B 0 = 0.1 T the minimum bubbling point is ε mb = 0.423 and U mb = 0.084 m/s, for B 0 = 0.15 T the mini-mum bubbling point is ε mb = 0.531 and U mb = 0.21 m/s. Notably, above B 0 = 0.188 T, the system is always stable since u d >u k at all flow rates.

CFD model
A set of test simulations were carried out to assess qualitatively the effect of the magnetic field on the fluidization regime. The aim is to verify the ability of the code in simulating whether the magnetic field is capable of suppressing the formation of bubbles in a bed that would normally be operating in the bubbling regime. In particular, a fluidized bed stabilized in LAST mode (i.e. magnetic field is applied to a freely bubbling fluidized bed) was simulated by CFD.
In the Eulerian-Eulerian approach to multiphase modeling, continuity and momentum balance equations in vector forms have to be solved for each phase : Of course, closure relations are also needed in order to properly model the particle phase and its interactions with the gas phase: for this purpose, the standard GKT model is adopted for estimating the rheological properties of the fluidized solid phase (Gidaspow, 1994) and standard drag models are adopted to estimate the momentum exchange between phases at the phase boundaries (as done by Brandani and Zhang, 2006). Complete details on the model equation implemented within the code can be found in the CFX documentation.
The additional forces due to the bed hydrodynamic elasticity (Eqns. 3-4: Brandani and Zhang, 2006;Busciglio et al. 2010) and the magnetic forces (Eqns. 14-15) were expressed as follows: The additional hydrodynamic forces, F f,y and F p,y , operate in vertical direction. Conversely, the magnetic body forces possess components in both the x and y directions. Maxwell's model equation has been used to obtain the additional terms reported in Eqns. 27-28. These are implemented in the code, since all the parameters except the void fraction gradients are constants.
As far as the numerical aspects are concerned, CFD simulations were performed in a 2D fashion, choosing a time step interval in the range between Δt = 10 -5 s and Δt = 10 -4 s and a computational grid consisting of 0.005 m square cells, with 200 cells along the vertical direction and 30 cells along the horizontal direction, thus resulting in a vertical extension of the domain equal to 1.0 m and a horizontal extension equal to 0.15 m. Grid independence was assessed, checking to make sure that cells only half the size (0.0025 m square cells) of those actually adopted did not change the simulation results.
The initial condition for the particle bed height is equal to 0.5 m, with the solid volume fraction within the bed set equal to 0.60. The lateral walls were modeled using the standard no-slip boundary condition. The upper section of the simulated geometry, or freeboard, was considered to be occupied only by gas. A simple pressure boundary condition was imposed at the top of the freeboard (i.e. fully developed flow condition). A Dirichlet boundary condition was employed at the bottom of the bed to specify uniform vertical gas inlet velocity throughout the distributor. Symmetry planes were imposed on the front and rear faces of the simulated bed in order to perform the simulation in a proper 2D fashion. Symmetry planes placed at the boundaries along the width direction of the computational domain causes all variables to be mathematically symmetric with thus no diffusion across the boundary except for the component of velocity normal to the boundary which is anti-symmetric. A typical Geldart type B system fluidized by air in magnetization LAST mode was simulated (d p = 300 μm, ρ p = 2500 kg/m 3 , χ p = 10). The fluid velocity was equal to 2.5 U mf , i.e. 0.22 m/s for the investigated system. The magnetic force is expected to cause particles to reduce their motion and form aggregates which are aligned along field lines. At high field intensity, the bed should eventually collapse, resulting in a fixed, anisotropic structure of particle aggregates.
A magnetic field, B, in the range between 0.05-0.6 T, was applied to a fully fluidized bed in order to stabilize the bed.
Each simulation was subdivided into the following three stages: • Starting from the initial steady condition of a settled bed, the gas inlet was simulated for 4s at Δt = 10 -4 s. During this time interval, the bed evolves from the initial condition to freely bubbling in pseudo steady state condition. Because of the inlet gas velocities adopted, a freely bubbling regime generally develops. • The freely bubbling condition at 4s is used as the initial condition for a second simulation stage, in which the magnetic field is present. Because of the very rapid transient dynamics expected, a reduced time step was adopted, i.e. Δt = 10 -5 s. In general, after 0.5 s subsequent to magnetic field application, the final steady state is practically achieved for the MSFB. • Finally, the simulation was carried out for a further 1s at a time step value of Δt = 10 -4 s: this allows simulation of the stabilized regime in a reasonable time. Typical CPU running times for this reference case were equal to about 4h for 1s of real time simulated with a fixed time step interval Δt = 10 -4 s on a Dell Dimension 8300 Personal Computer.

Data analysis
Data obtained from the simulations were post-processed in order to quantify the effect of the superimposed magnetic field. The following instantaneous parameters were chosen: • the instantaneous polydispersity index of volume fraction distribution within the bed, PI, was chosen to characterize the uniformity of bed voidage. This parameter is defined as: where the general moment ratio ϕ ij (t) is defined on the basis of the solid volume fraction values of the com-putational cells placed below the initial bed height H 0 (t = 0). Notably, it was decided to analyze all the computational cells below the initial bed height to avoid measurement of the properties of the freeboard region. The larger the polydispersity index PI(t), the less uniform is the instantaneous volume fraction distribution inside the bed. • the instantaneous mean value of bed height, H 0 , was chosen to characterize bed expansion below the initial bed height. The freeboard extent is computed as the average height if the upper domain region has a solid volume fraction below 0.15 (this is equivalent to the conventional transition value between the emulsion phase and bubble phase in fluidized beds).

Results
The first qualitative evaluation of the simulation ability to capture the typical behavior of a MSFB can be performed easily by observation of the particle volume fraction map sequences reported in Fig. 2. The maps report data about the glass ballotini simulation, d p = 300 µm; fluidized at U = 2.5 U mf with a magnetic field intensity ranging from B = 0.05 T to B = 0.6 T.
Taking as an example the sequence simulated at B = 0.2 T (Fig. 2b), which is immediately above the minimum stabilization field, it is possible to observe that the sequence embraces all stages from the initial condition to fully developed bubbling regime conditions (from 0.0 to 4.0 s), and the subsequent transitory due to application of the magnetic field (i.e. from 4.0 to 5.5 s) until steady state conditions are achieved.
It is worth noting that the code is able to simulate bubble formation correctly, at least in a qualitative way, and growth along the bed appears to be correctly predicted. Notably, after application of the magnetic field, the fluid dynamic regime suddenly changes toward a stabilized state. Under these conditions, all bubbles disappear in all cases simulated.
As it can be observed in Fig. 2a, when the magnetic field is set at values lower than 0.2 T, the stabilizing effect is able to partially damp out bubbles within the bed, but local in homogeneities are still evident. In this case, the applied field intensity was not sufficient to achieve a fully stabilized state. In all likelihood, the transition from partial to full bed stabilization can be found at magnetic field values close to B = 0.2 T as discussed earlier. In the case of higher field values (B = 0.4 and 0.6 T shown in Fig.  2d-e), the bed is fully stabilized but a longitudinal voidage gradient appears parallel to the field direction. Notably, this increases the bed volume occupied by the solid phase and the magnetic field also increases.
This effect can be theoretically predicted taking into ac-count the one-dimensional momentum balances (Eqns. 1-2). At equilibrium, if there are no gradients in the void fraction of the bed, the bed weight must equate the drag force, hence it is possible to calculate the void fraction of the homogeneous bed, which for U = 0.22 m/s yields ε = 0.537. Most of the bed will be at this condition, except in the vicinity of the freeboard where a gradient in the void fraction is inevitable. At equilibrium, we can assume that u p = 0, u f = U/ε and all time gradients are nil. If the fluid and particle momentum balance are combined to eliminate pressure, the following differential equation is obtained: Eqn. 30 clearly shows that there will be a region between ε 0 = 0.537 and ε 1 = 1-where the RHS of the equation is zero-where there will be a transition. Eqn. 30 can be used to estimate the stratification at the freeboard as a function of the magnetic field as shown in the figure below (Fig. 3) Clearly, further computational and experimental investigations are needed to investigate this behavior in detail, but the full simulations are clearly in qualitative agreement with the predicted effect. In order to quantify and compare the simulation results after the change in the applied magnetic field, in Fig. 4a the time evolution of the polydispersity index is reported for all the simulations shown in Fig. 2.
Before the onset of the magnetic field (t = 0-4 s), it is possible to observe high PI values with a markedly oscillating behavior due to the bubbles rising through the bed. After application of the magnetic field, a sudden decrease of the polydispersity index underlines the stabilization of the system for all cases investigated.
Notably, the effect of the applied magnetic field is practically instantaneous. Moreover, the quantitative analysis of PI shows the stabilizing effect that the magnetic field has on the bubbling bed: in the case of B < 0.2 T, the PI value decreases but still oscillates after the magnetic field onset, while higher magnetic field intensities give rise to a smooth and stable PI value. This is further confirmed by Fig. 4b, in which time-averaged PI values are reported as a function of the superimposed magnetic field intensity. It is possible to observe that at B values below 0.2 T, high polydispersity values are found because of incomplete bubble suppression, while at values higher than B = 0.4 T, even if complete bubble suppression occurred, the PI value increases slightly because of bed stratification. Hence, these simulations confirm the ability of the model to predict the existence of a critical B value to achieve a fully stabilized bed.
Clearly, the onset of a diffused freeboard (qualitatively observed in Fig. 2 and further discussed in Fig. 3) can be quantified by showing the bed height H 0 as a function of time (as reported in Fig. 5a) and relevant time-averaged values (reported in Fig. 5b). Notably, the increase in the average bed height becomes important at a high magnetic field intensity (especially above B = 0.4 T).
Concluding, it is possible to point out that even if some limits are still present, because of intrinsic limitations of the multifluid approach (Eulerian-Eulerian) and magnetic model simplification, the CFD model adopted is able to capture the main characteristics of the MSFB with reasonable accuracy and within acceptable simulation times from an engineering point of view.

Conclusions and future aims
An Eulerian-Eulerian CFD model to simulate the behavior of the MSFB was successfully developed. The CFD model correctly predicts the onset of different stabilization levels depending on the intensity of the magnetic field applied and provides further insights on the behavior of MSFBs.
Even if work still has to be done with regard to the implementation of more comprehensive mathematical models (including inter-particle forces as an example), this result may be regarded as very promising, being among the first CFD simulations of MSFBs with Eulerian-Eulerian models.