CFD Modelling of Coupled Multiphysics-Multiscale Engineering Cases

Abstract

Many of the engineering problems have multiphysics and multiscale nature. Non-isothermal flows, stirred reactors, turbulent mixing and membrane filtration, are prevalent cases in which the coupling of several physics phenomena is required for the adequate prediction of overall behaviors. Also, a multiscale analysis, where the same phenomenon is analyzed at different scales, can lead to better understanding of the phenomena, which can be used in optimization and to provide adequate scale-up methodologies. Studies incorporating both multiscale and multiphysics analysis are rarely addressed in literature; in fact, these kinds of problems will be the research challenge in the next years. Computer fluid dynamics (CFD) techniques have shown to be promising to deal with these kinds of systems. In this chapter, these are used to implement a multiscale analysis of the hydrodesulphurization (HDS) process for light gas-oil (LGO). The aforementioned is carried out by the analysis of mass an energy transport at: (1) microporous (MP) scale, (2) pseudo-homogeneous catalyst (PHC) scale, and by analysis of (3) momentum and mass transport at reactor scale (RS). In addition, a particular discussion is made regarding the proper establishment of the model, its validation, the use of different boundary conditions, its justification; and the dependence of solutions of parameters and initial and boundary conditions.

1. Introduction

The multiscale phenomena are common issues for many applications in several processes and phenomena of interest overall engineering fields. From the design, implementation and optimization of processes and equipment, to the Process Safety Engineering, this multiscale nature represents a great challenge to overcome. An example of this is the study of industrial fires, in which the fire can reach over 50 m in height, and the produced smokes may have effects on the environment in lengths of tens or even hundredths of kilometers. Also, in these phenomena, the hydrodynamics of circulating air has a crucial effect in the behavior of the fires, which have a different characteristic length scale, of a few kilometers. Furthermore, it is known that there is also a turbulent transport of energy, which has a predominant contribution over the energy transport in these phenomena, and important parameters as the turbulent kinetic energy varies in the inertial subrange of 1 × 10⁻³ m; while the scale for the dissipation of kinetic energy (Kolmogorov scale) is about 1 × 10⁻⁴ m. An even smaller scale occurs in the molecular diffusion of the species involved in the phenomenon (Batchelor scale), which is in the range of 1 × 10⁻⁵ − 1 × 10⁻⁶ m.

Another example of this multiscale nature is present in multiphase catalytic chemical reactors, where three different length-scale analysis levels can be distinguished: (1) reactor level, (2) catalyst or catalyst’s clusters level, and (3) catalyst microstructure level. Every level prior mentioned having physical and chemical phenomena taking place at different scales, and being of individual interest for different kinds of analysis. For example, at the last-mentioned level, there is an interest into investigating over the size and distribution of noble metal crystals, commonly named adsorbed islands, their interaction with the supporting material, and other key aspects that determine the catalyst activity. These phenomena take place at scales around 1 × 10⁻⁹ − 1 × 10⁻⁸ m; it is noteworthy that in the study of this low scale, even the continuum assumptions are not valid, and therefore there is a need to implement molecular analysis. At the catalyst cluster lever, the characteristic length can be set to the pellet radius, which are found commercially in the range of 1 × 10⁻³ − 1 × 10⁻² m. In this scale, there is an interest in the phenomena taking place inside the porous microstructure of the catalyst and in the boundary layer, where coupled phenomena of mass and energy transport occur in lengths below 1 × 10⁻³ m, and phenomena of resistance to the transport between the phases occur below 1 × 10⁻⁴ m. In the large scale, the reactor length, the phenomena occur in length scales of around 1 × 10² m. In this, there is an interest in the study of maldistribution, incomplete wetting, phases distribution, hot spot formation, among other issues, which can be caused by phenomena at lower scales. For both of the mentioned examples, a similar discussion can be made regarding the time scales. Figure 1 depicts the length scale of the phenomena taking place in the study of catalytic reactor. This example will be further developed in order to give insight into the multiscale and multiphysic nature of engineering problems.

Figure 1.

Schematic representation of the multiscale in catalytic reactors.

In addition, it is important to note that both examples also present a strong multiphysic nature; both involve chemical reactions of hundreds, up to thousands, of chemical species, multiphase transport of species, energy transport generated by the reaction nature, and momentum transport of fluid phases that interact with solids. Due to these complexities, mathematical methods have been developed in order to better understand these interactions and the phenomena that is hardly observable by experimental methods. In this context, CFD techniques have come into view as promising alternatives into the multiphysic-multiscale modelling. However, these approaches face the obstacle of the high nonlinearity and high grade of coupling between the phenomena, physics and scales. Furthermore, in the catalytic reactor study, further nonlinearities are present as a result of the nonlinear dependence of the concentrations and temperature fields to the chemical reactions. Additionally, in the reactor scale, more coupling nonlinearities are present [1] due to the extra phenomena, such as adsorption/desorption at solid-fluid interfaces, caused by heterogeneities of the porous media. Then, it is natural to think that a precise description of the effects and dependence between the scales must be taken into account to properly describe the phenomena.

In scientific community, it has been recognized that in complex systems, the dominant issue is the existence of multi spatial-temporal scales and its multiphysics nature; and it has identified that its comprehension is at the frontier of state of art of process engineering, as well as in science and technology of many fields and disciplines [2]. In spite of the above, the knowledge developed in engineering fields until the first years of the present century has been extending from understanding macro-scale field (mass, energy, velocity, etc.) distribution and individual phenomena, focusing on mechanisms at microscale, to the understanding of coupling between phenomena at different scales. To enable these studies, it has been necessary to make average assumptions over fields transport equations to studying heterogeneity both in time and in space. However, due to the limitation in the knowledge to deal with non-equilibrium and nonlinear phenomena, the quantification of chemical and physical processes at different scales is yet a major challenge. The use of average and linear approaches and simplification of heterogeneities is insufficient to deal with multiscale spatial-temporal structures and nonlinear phenomena, then the multiscale analysis have emerged as promising tool to improve models, theories, and knowledge about these systems. The development of computers and advances in numerical methods have allowed the development of software and advances in measurement at small and non-invasive scales, and have also favored the development of the multiscale analysis [3, 4, 5, 6, 7].

A critical step in the development of multiscale models is the establishment of bridges that allow the coupling of variables and parameters between the different scales, and the transfer of information between these scales [4]. Regarding this latter aspect, it is also important to determinate if the transferred information between the scales is required to be in two-ways (up to down, and down to up scales) or one-way (down to up scales) [8]. The consideration of these aspects is what will allow the analysis of complex systems moving from the reductionism of average simplifications to a holistic form of study.

Great efforts into the CFD analysis of multiscale systems have been developed in recent years. It is noteworthy the work of Ding and co-workers [7], who developed a multiscale methodology aided by CFD simulation, for a catalytic distillation with bale packings, considering micro and macroscales, focusing only in the hydrodynamics. In the microscale, a volume of fluids (VOF) method is implemented in a representative elementary unit (REU) to simulate the gas-liquid flow considering the packing geometry, in order to study the liquid split proportion. While in the macroscale, a unit network model is developed, where the liquid split proportion from the microscale model is taken as an input, in order to measure the liquid holdup and pressure drops. It is important to note that in this work, the two scales are communicated in only one way. This means that the information of the microscales phenomena scales up to the macroscale, but the effects from the macroscale to the microscales are neglected. Also, no multiphysic analysis in developed. Several works where multiscale analysis is carried out only in the hydrodynamics can be found in literature [9, 10]; however, works where coupled multiscale and multiphysic analyses are performed, are scarce.

As an example of the aforementioned, the works of Xie and Luo can be mentioned [11]. They developed an Euler-Euler two-phase model and a population balance model (PBM) to simulate a liquid-liquid suspension polymerization process. The model considers momentum, turbulence, mass and energy transport equations, while PBM considers a Breakage Kernel and Coalescence Kernel. It is important to highlight is that these differential equations have the same domain, implying solutions at the same time and length scales. The authors consider that as the polymerization process occurs in a lower length-scale, and the kinetics of this process occurs in an atomic scale, the multiscale is captured through the effects of these. These kinds of models have been recently developed, and there is still a need to clarify certain aspects regarding the multiscale-multiphysic nature, and how should it be addressed and incorporated to the models. Thus, further works that incorporate the multiscale-multiphysic analyses, and address the establishment of bridges that connect the scales and how these communicate the information, are yet desirable.

In this work, an effort to contribute to the knowledge of multiphysic-multiscale modelling is developed. The analysis of a hydrodesulphurization (HDS) process for a light gas-oil (LGO) through COMSOL Multiphysics CFD simulations is addressed and discussed, considering three different scales (i) micropores (MP) scale, (ii) pseudo-homogeneous catalyst (PHC) scale, and (iii) reactor scale (RS). The analysis takes into account the effect of the microstructure geometry on upper scales, and how is this information captured and communicated. In addition, the effect of the reactor scale over the lower scales is analyzed, by means of the observed differences when the catalyst model is solved without the effect of the reactor behavior. By the last, a discussion regarding the boundary conditions establishment, model validation, and dependence of solution of the model parameters, is presented.

2. HDS process and HDS reactor

The hydrodesulphurization (HDS) is part of catalytic hydrotreatment (HDT) process where the content of some crude oil contaminants containing sulfur is reduced using hydrogen over a catalyst of NiMo or CoMo supported on Al₂O₃. This is one of the most important processes in crude oil refining, because it allows reducing the emission of SO_x and NO_x, which are synthesized by fuel combustion. These emissions are strong environmental contaminants; in addition, they can prejudice the performance of the catalysts used in refining processes, as well as the catalysts used in catalytic converters of vehicles [12].

The most important equipment in a HDT unit is the three-phase HDS reactor, and thus this has been target of extensive investigation and modelling. In HDS reactors, gas and liquid phases (hydrogen and a fraction of hydrocarbons) are contacted with a solid phase (catalyst). The reactions occur between the dissolved gas reactant and the liquid-phase reactant at the surface of the catalyst. Depending on whether the main mass-transfer resistance is located, three-phase catalytic fixed-bed reactors can operate either, with a continuous gas and a distributed liquid phase (trickle operation), or with a distributed gas and a continuous liquid phase (bubble operation). Commercial HDS reactor usually operates in a trickle-bed regime, with concurrent downward flow of gas and liquid over a randomly fixed bed of catalyst particles where the reactions take place [13]. The HDS reactor operates at elevated pressures and temperatures because the solubility of gases in liquids increases with rising partial pressure and the reaction rate is favored with those temperature magnitudes.

Given the importance of the HDS process and the growing pressure exerted by the new environmental regulations reducing the maximum limits of sulfur content in fossil fuels, have led to the need for optimization of HDS reactors; which in turn requires deeper and more detailed knowledge of the phenomena that occur in the HDS reactor. To obtain the knowledge that allows the optimization of the HDS reactor, the CFD techniques are very promising, as they allow to obtain punctual and average values for the fields of concentration, temperature, etc.

In reactor modelling, a common approach is to make corrections to the intrinsic reaction rate (r_i) to take into account the effects of the mass and energy transport resistances inside and outside the catalyst, through a called effectiveness factor (η), defined by Eqs. (1) and (2).

${\left(r_i\right)}_{\mathit{\text{observed}}}=\eta {\left(r_i\right)}_{\mathit{\text{intrinsic}}}$E1

$\eta =\frac{{\int }_{V_{\mathit{\text{catalyst}}}}{r}_i\left({c_i}^{\text{'}}s,T\right)\mathit{dV}}{{\int }_{V_{\mathit{\text{catalyst}}}}{r}_i\left({c_i}^{\text{'}}s,T\right){|}_{\mathit{at surface condition}}\mathit{dV}}$E2

The correction is due to that “it is not possible or it is hard” to solve the exact heterogeneous mass and energy transport equations with superficial reaction, which are valid inside the pores and the solid matrix of catalyst; the aforementioned, is due to that the geometry of the equations domains is unknown or very complicated. Therefore, the mass and temperature fields necessary to evaluate the intrinsic reaction rate in Eq. (1), are evaluated through solving pseudo-homogeneous mass and energy transport equations with reaction taking place in all catalyst. In this pseudo-homogeneous model, the porous nature and its mass and energy transport characteristics, inherent of the real heterogeneous catalyst, are represented or incorporated through two effective transport coefficients, the effective diffusivity and conductivity coefficients (D_eff, K_eff). It is important to emphasize that in the described approach, there are two events where the transport of information from a small to a larger scale is carried out. In precise form, the first take place when the effective transport coefficient captures geometry characteristics and transport phenomena happening at porous microstructure scale; and the second takes place when mass and energy features, inside catalytic particle are incorporated to correct the reaction rate, which is used at reactor scale.

In this chapter, an example where the multiphysics and multiscale nature of trickle bed reactor’s (TBR’s) for HDS process are discussed is addressed. The CFD models are constituted by models at three different scales: (a) porous microstructure scale of catalyst, (b) catalyst particle scale, and (c) reactor scale. In Figure 1, the three scales and the geometric details of the implementing models here are shown.

3. Multiscale CFD models

3.1. Micropores CFD model

At MP scale, the heterogeneous mass and energy transport equations are solved, taking into account the transport phenomena occurring inside pores and solid matrix of a catalytic particle. In this system, the effective transport coefficients are evaluated.

3.1.1. Porous microstructure geometric model

To perform the mentioned above, a representative model of the geometry microstructure of a catalyst was constructed using a vectored model of a real porous media taken from a micrograph found in literature [14]. This was used and adapted in order to replicate parameters comparable with typical values of HDS systems such as pore diameter d_p = 20 − 200 nm [15] and porosity 0.3 ≤ ε_λ ≤ 0.6 [16]. It is worth noticing that a real geometric model of pores distribution has been built for the catalytic particle of 0.35 mm of diameter, at scales of 4:1 and 2:1 of the real scale.

Figure 2 shows detail of geometric representation of catalytic particle incorporating an explicit porous microstructure, which is constituted by a solid matrix (phase-σ) and the fluid phase formed by the interstitial spaces left by the solid matrix (phase-λ). In addition, the solid-fluid interphase (A_σλ) where the superficial reaction takes place, is illustrated.

Figure 2.

Representation of the heterogeneous catalyst micropores model.

3.1.2. Heterogeneous mass and energy transport model

The transport model for the MP model, is formed by local mass and energy transport equations. For the case of the mass transport, it is considered that there is only transport in the interstitial phase, because the porous matrix is considered to be impermeable; while for energy transport both phases carry energy, so there is an equation for each domain.

It is important to note that due to capillarity forces, the micropores are completely filled [12], consequently the fluid inside of catalyst’s pores is stagnant, there are no convective contributions to the mass and energy transport, and there is no need for a momentum balance equation. Therefore, the mass and energy transport models are reduced to Eqs. (3) and (4), respectively [17]

$\nabla \cdot \left({\mathrm{D}}_i^{\lambda }\nabla C_i^{\lambda }\right)=0;i=\left\{R-S,H_2,H_{2}S\right\}\mathit{\text{specie}}$E3

$-k_j{\nabla }^2{T}^j=0;j=\left\{\sigma ,\lambda \right\}\mathit{\text{phase}}$E4

It is important to note that in the transport Eqs. (3) and (4) there is no reaction term and heat source term, respectively, thus concentrations and temperature fields are due to the non-homogeneous boundary conditions. In fact, the superficial reaction, as well as the heat generation by reaction are considered as a no-homogeneous boundary condition (see Table 1).

Mass transfer (i ‐ specie = R ‐ S, H2, H2S)
$-{\mathbf{n}}_{\mathit{\sigma \lambda }}\cdot {\mathrm{D}}_i^{\lambda }\text{∇}{C}_i^{\lambda }={\nu }_i{a}_v{r}_{\mathit{HDS}}$ (superficial reaction)	at Aσλ	(5)
$C_i^{\lambda }={\left.C_i\right|}_{\mathit{\text{bulk}}}$ (bulk phase concentration)	at Aλe	(6)
Heat transfer in liquid phase (phase-λ)
−nσλ ⋅ kλ∇Tλ = av(−ΔHHDS)rHDS (energy generation at porous interphase)	at Aσλ	(7)
Tλ = T0 (bulk phase temperature)	at Aλe	(8)
Heat transfer in solid phase (phase-σ)
−nσλ ⋅ kσ∇Tσ = av(−ΔHHDS)rHDS (energy generation at porous interphase)	at Aσλ	(9)
Tσ = T0 (bulk phase temperature)	at Aσe	(10)

Table 1.

Boundary conditions of heterogeneous transport CFD model.

When it is considered that the reaction takes place in all the fluid-solid interface within the pores of the catalyst, a simplification has been made because actual the reaction takes place in a lower scale on the active phase of the catalyst, which are nothing more than small cumulus of crystals of metals or metal salts, which are deposited on the surface (so-called adsorbed islands). Moreover, on these adsorbed islands phenomena of adsorption and desorption of products and reagents take place.

The foregoing means that both reaction and adsorption/desorption phenomena occur at discrete locations on the surface of the catalyst and not on the entire surface. Fortunately, the experimentally determined kinetics expressions takes into account the proper considerations in order to express the reaction that takes place in punctual as a reaction occurring throughout the catalytic interfacial surface. Also, the predictions of the models considering surface reaction have proven to be accurate enough [18].

The boundary value problem specified by Eqs. (3) and (4) are set to satisfy the boundary conditions shown in Table 1.

Where A_λe and A_σe represent the external catalyst boundaries; −n_σλ is the unitary normal vector that points from σ ‐ phase to λ ‐ phase; ν_i is the stoichiometric coefficient for each specie; and C_i|_bulk is the concentration of the specie i at reactor bulk conditions.

Note that in the boundary conditions (5), (7), and (9), the reaction rate (r_HDS) is multiplied by the parameter a_v, which is defined as the ratio of interstitial volume to fluid-catalyst interfacial area (a_v = V/A_σλ) [18]. This parameter appears by virtue of averaging process of mass transport equations and superficial reaction rate and allows us to relate the reaction rate obtained by experimental data ${⟨r_i^{\omega }⟩}^{\omega }\left[\mathrm{mol}/{\mathrm{m}}^3\mathrm{s}\right]$ with the superficial velocity r_i[mol/m²s].

$r_i=a_v{⟨r_i^{\omega }⟩}^{\omega }$E11

3.2. Effective coefficients evaluation

A theoretical development to evaluate these coefficients has been established by Whitaker and co-workers, within the framework of the method of volume averaging [18]. The expressions that allow the evaluation of the effective diffusivity and effective conductivity coefficients are Eqs. (12) and (13), respectively. Note that these expressions are valid for any engineering case involving diffusion and conductive heat transfer through porous media where reaction takes place at interphase area.

${\mathbf{D}}_{\mathit{eff},i}={\mathrm{D}}_i^{\lambda }\left(\mathbf{I}+\frac{1}{V_{\lambda }}\underset{{\mathrm{A}}_{\mathit{\lambda \sigma }}}{\int }{\mathbf{n}}_{\mathit{\sigma \lambda }}{\mathbf{b}}_{\lambda }\mathit{dA}\right)$E12

$\frac{{\mathbf{K}}_{\mathit{eff}}}{k_{\lambda }}=\left({\epsilon }_{\lambda }+{\epsilon }_{\sigma }\kappa \right)\mathbf{I}+\frac{\left(1-\kappa \right)}{\mathtt{V}}\underset{{\mathrm{A}}_{\mathit{\lambda \sigma }}}{\int }{\mathbf{n}}_{\mathit{\lambda \sigma }}{\mathbf{b}}_{\Lambda }\mathit{dA}$E13

In the equations above, ε_i represent the volume fraction occupied by each i-phase; I is the identity tensor; and κ is the quotient of conductivities of the solid phase to fluid phase κ = k_σ/k_λ.

To evaluate the effective coefficients, the evaluation of so-called closure vectors b_λ and b_Λ are required, which is done through the solution of boundary value problems described by the following equations.

${\nabla }^2{\mathbf{b}}_{\lambda }=0$E14

${\nabla }^2{\mathbf{b}}_i=0;i-\mathit{\text{phase}}=\Lambda ,\sigma$E15

These boundary value problems are the result of the averaging process for the punctual mass and energy transport Eqs. (3) and (4), of a decomposition of scales, and of a proposal for a solution of the deviations field problem using the source terms boundary value problem for deviations for mass and energy average equation. For more details, the reader is invited to review the work of Whittaker [18] and extensive literature regarding Method of volume averaging [19, 20, 21].

On the other hand, in Eqs. (14) and (15), the subscripts Λ and λ, both refer to the fluid phase, however, it is important to differentiate them as they come from the solution of a different boundary value problem.

The boundary conditions set to both problems are shown in Table 2.

Mass transfer
−nλσ ⋅ ∇bλ = nλσ	at Aσλ	(16)
bλ(r + li) = bλ(r)	for i = 1, 2, 3, …	(17)
Heat transfer
bΛ = bσ	at Aσλ	(18)
−nλσ ⋅ ∇bΛ =  − nΛσ ⋅ κ∇bσ + nΛσ(1 − κ)	at Aσλ	(19)
bΛ(r + li) = bΛ(r);   i − phase = Λ, σ	for i = 1, 2, 3, …	(20)

Table 2.

Boundary conditions for closure boundary value problems.

In these, r is the position vector that locates any points in the average volume, and l_i represent the three non-unique lattice vectors that are required to describe a spatially periodic porous medium [22]. Then, boundary conditions (17) and (20) are actually periodicity conditions, implying that these boundary value problems are usually solved in periodic domains inherent to periodic representative unitary cells (RUC).

It must be noted that the boundary value problem for the closure vectors are essentially geometrical, and that the vector field is generated by non-homogeneous boundary conditions; this means that the vector field is generated by the presence of solid-fluid interphases (A_σλ) inside of porous media. Also, it should be noted that information of the porous structure and its effects are captured by the closure vectors, and that the closure vectors for mass and energy are different. Implying that something else that the geometrical characteristics are captured through them. This information passes to upper scales through the effective transport coefficients (Eqs. (12) and (13)), implying that the porous structure affects the mass and energy transport at pellet scale. Then, in can be seen that this development is setting a bridge between the different length-scales.

3.3. Pseudo-homogeneous CFD model

The mass and energy pseudo-homogeneous transport equations for a catalytic particle are well established in the literature [23], for the analyzed HDS reaction here, these take the form:

$\nabla \cdot \left({\epsilon }_{\lambda }{\mathbf{D}}_{\mathit{eff}}\nabla {⟨C_i^{\omega }⟩}^{\omega }\right)={\nu }_i{r}_{\mathit{HDS}}i-\mathit{\text{species}}=R-S,H_2,H_{2}S$E21

$\nabla \cdot \left({\mathbf{K}}_{\mathit{eff}}\cdot \nabla {⟨T⟩}^{\omega }\right)=\left(-\mathit{\Delta H}\right)r_{\mathit{HDS}}$E22

In these expressions, ${⟨C_i^{\omega }⟩}^{\omega }$ and ⟨T⟩^ω represent the intrinsic average of concentration for the i-specie and temperature fields, which are quite different from the punctual concentration and temperature fields $\left(C_i^{\lambda },T^j\right)$ of the heterogeneous model described by Eqs. (3) and (4). To clarify the above it is important to remember that the average fields of concentration and temperature present changes in lengths of scale of order of l_pellet < D_pellet, whereas the point fields undergo changes in the length scale of the order of l_porous < d_pore.

On the other hand, the effective coefficients of transport (D_eff, i, K_eff) necessary to solve the model at catalyst scale can be obtained by experimental data or by theoretical approaches. In the present case, these were obtained by the solution of boundary value problem for closure vector (Eqs. (14)–(20)) in conjunction with Eqs. (12) and (13).

The appropriate boundary conditions for the transport of matter and energy in the catalyst are shown in Table 3.

Mass transfer (i − specie = R − S, H2, H2S)
$\text{∇}{⟨C_i^{\omega }⟩}^{\omega }=0$ (concentration continuity)	at r = 0	(23)
${⟨C_i^{\omega }⟩}^{\omega }={\left.C_i\right|}_B$ (bulk phase concentration)	at r = Rp	(24)
Heat transfer
∇⟨T⟩ = 0 (temperature continuity)	at r = 0	(25)
⟨T⟩ = ⟨T⟩B (bulk phase temperature)	at r = Rp	(26)

Table 3.

Boundary conditions for pseudo-homogeneous transport model.

As can be observed, temperature and concentration of each specie at the bulk fluid inside reactor are required; this implies that there is a need to transfer information between the phenomena that take place on the catalytic scale and on the reactor scale. The transfer of the information of the nature of the resistances on the transport at catalyst scale, toward reactor scale, can be carried out through the effectiveness factor, as described in Section 2.1.

Figure 3 shows the relationship and differences between heterogeneous and pseudo-homogeneous mass transport at catalytic particle.

Figure 3.

Concentration field for a catalytic particle obtained through the transport models at porous and catalyst scales.

3.4. Heterogeneous CFD reactor scale model

The CFD reactor model considers an Eulerian approach, where gas and liquid phases are considered as interpenetrating, implying that both fluid phases have been the same domain. In the case of the solid phase, an explicit geometry for a fixed bed was built considering spherical catalytic particles of 0.35 mm of diameter. Figure 4 shows the geometrical details for both interstitial (γ, β ‐ fluid phases) and catalytic (ω ‐ phase) domains.

Figure 4.

Geometrical details of the TBR model.

The complete CFD model for the three-phase reactor consists of two momentum transport equations, one for each fluid phase, eight mass transport equations, one per specie, and three closures for the interaction between the solid, gas, and liquid phases.

3.4.1. Hydrodynamic model for TBR

The hydrodynamic model at steady state is constituted by continuity and momentum transport equations:

${\epsilon }_{\gamma }{\rho }_{\gamma }\nabla \cdot {\mathbf{v}}_{\gamma }=0\left({\rho }_{\gamma }=\mathit{cte}\right)$E27

${\epsilon }_{\beta }\nabla \cdot \left({\rho }_{\beta }{\mathbf{v}}_{\beta }\right)=0$E28

${\epsilon }_{\gamma }{\rho }_{\gamma }\left({\mathbf{v}}_{\gamma }\cdot \nabla \right){\mathbf{v}}_{\gamma }={\epsilon }_{\gamma }\nabla \cdot \left[-P\mathbf{I}+{\mu }_{\gamma }\left(\nabla {\mathbf{v}}_{\gamma }+{\left(\nabla {\mathbf{v}}_{\gamma }\right)}^T\right)\right]+{\mathbf{F}}_{\gamma }/{\epsilon }_{\gamma }+{\rho }_{\gamma }\mathbf{g}$E29

${\epsilon }_{\beta }{\rho }_{\beta }\left({\mathbf{v}}_{\beta }\cdot \nabla \right){\mathbf{v}}_{\beta }={\epsilon }_{\beta }\nabla \cdot \left[-P\mathbf{I}+{\mu }_{\beta }\left(\nabla {\mathbf{v}}_{\beta }+{\left(\nabla {\mathbf{v}}_{\beta }\right)}^T\right)-\frac{2}{3}{\mu }_{\beta }\left(\nabla \cdot {\mathbf{v}}_{\beta }\right)\mathbf{I}\right]+{\mathbf{F}}_{\beta }/{\epsilon }_{\beta }+{\rho }_{\beta }\mathbf{g}$E30

In these, v_i and ε_i are the local interstitial velocity and volume fraction for both gas and liquid phase, respectively. In Eqs. (27) and (28), the gas phase is considering as compressible fluid.

The term (F_i/ε_i) takes into account the momentum exchange between the phases through the momentum exchange coefficient K_ij, which has three contributions: liquid-gas, liquid-solid, and gas-solid interactions.

${\mathbf{F}}_i=\sum _{j=1}^n{K}_{\mathit{ji}}\left({\mathbf{v}}_j-{\mathbf{v}}_i\right)$E31

In the CFD model, the Attou momentum exchange model, which is actually a closure for the momentum transfer equations, was incorporated [24], which is considered adequate to take account of the interaction between phases in fixed bed three-phase reactors.

$K_{\mathit{\beta \gamma }}={\epsilon }_{\beta }\left[\frac{E_1{\mu }_{\beta }{\left(1-{\epsilon }_{\beta }\right)}^2}{{\epsilon }_{\beta }^2{d}_p^2}{\left(\frac{{\epsilon }_{\omega }}{1-{\epsilon }_{\beta }}\right)}^{0.667}+\frac{E_2{\rho }_{\beta }\left({\mathbf{v}}_{\beta }-{\mathbf{v}}_{\gamma }\right)\left(1-{\epsilon }_{\beta }\right)}{{\epsilon }_{\beta }{d}_p}{\left(\frac{{\epsilon }_{\omega }}{1-{\epsilon }_{\beta }}\right)}^{0.333}\right]$E32

$K_{\mathit{\beta \omega }}={\epsilon }_{\beta }\left[\frac{E_1{\mu }_{\beta }{\left(1-{\epsilon }_{\beta }\right)}^2}{{\epsilon }_{\beta }^2{d}_p^2}{\left(\frac{{\epsilon }_{\omega }}{1-{\epsilon }_{\beta }}\right)}^{0.667}+\frac{E_2{\rho }_{\beta }{\mathbf{v}}_{\beta }\left(1-{\epsilon }_{\beta }\right)}{{\epsilon }_{\beta }{d}_p}{\left(\frac{{\epsilon }_{\omega }}{1-{\epsilon }_{\beta }}\right)}^{0.333}\right]$E33

$K_{\mathit{\gamma \omega }}={\epsilon }_{\gamma }\left[\frac{E_1{\mu }_{\gamma }{\epsilon }_{\omega }^2}{{\epsilon }_{\gamma }^2{d}_p^2}+\frac{E_2{\rho }_{\gamma }{\mathbf{v}}_{\gamma }{\epsilon }_{\omega }}{{\epsilon }_{\gamma }{d}_p}\right]$E34

where E₁ and E₂ are de Ergun constants, and μ_i and ρ_i are i-phase viscosity and density, respectively.

3.4.2. Mass transport model for i-specie at TBR

The three-phase mass transport model is constituted by the following set of differential equations:

$\nabla \cdot \left(-{\mathrm{D}}_i^{\beta }\nabla ⟨C_i^{\beta }⟩\right)+{\mathbf{v}}_{\beta }\cdot \nabla ⟨C_i^{\beta }⟩=N_i^{\mathit{\beta \gamma }}/{\epsilon }_{\beta }i-\mathit{\text{specie}}=H_2,H_{2}S$E35

$\nabla \cdot \left(-{\mathrm{D}}_i^{\gamma }\nabla ⟨C_i^{\gamma }⟩\right)+{\mathbf{v}}_{\gamma }\cdot \nabla ⟨C_i^{\gamma }⟩=N_i^{\mathit{\gamma \beta }}/{\epsilon }_{\gamma }i-\mathit{\text{specie}}=R-S,H_2,H_{2}S$E36

$\nabla \cdot \left({\epsilon }_{\lambda }{\mathrm{D}}_{\mathit{eff},i}\nabla {⟨C_i^{\lambda }⟩}^{\lambda }\right)={\nu }_i{r}_{\mathit{HDS}}i-\mathit{\text{specie}}=R-S,H_2,H_{2}S$E37

For case of fluid phases, the mass transport considers both the convective and diffusive contribution, and the chemical reaction does not take place in these phases. In addition, both fluid phases are coupled by a volumetric flux exchange term $\left(N_i^{\mathit{\beta \gamma }}/{\epsilon }_{\beta }\right)$ defined by Eq. (38), in which $K_{g,i}^{\mathit{\beta \gamma }}$ is the gas-liquid mass transfer coefficient.

${\mathbf{N}}_i^{\mathit{\beta \gamma }}=K_{g,i}^{\mathit{\beta \gamma }}\left(⟨C_i^{\beta }⟩R_{g}T/H_i-⟨C_i^{\gamma }⟩\right)$E38

For the catalytic phase, volumetric chemical reactions, and only the diffusive contribution to the mass transport are considered.

Table 4 shows the boundary conditions adequate for mass and momentum transport for the TBR model for the HDS process.

Hydrodynamics i − phase = β, γ
${\mathbf{v}}_i=-\mathbf{n}{v}_i^0$ (inlet velocity)	at z = LR	(39)
P =  − P^0; n[μi(∇vi + (∇vi)^T)] = 0 (oulet pressure)	at z = 0	(40)
vi = 0 (no slip condition)	at r = rR and Aωi	(41)
$\mathbf{n}\cdot \mathbf{v}=0;{\mathtt{K}}_i-\left({\mathtt{K}}_i\cdot \mathbf{n}\right)\mathbf{n}=0$ (symmetry)	at plane in x = 0	(42)
${\mathtt{K}}_{\gamma }={\mu }_{\gamma }\left(\nabla {\mathbf{v}}_{\gamma }+{\left(\nabla {\mathbf{v}}_{\gamma }\right)}^T\right)\mathbf{n};{\mathtt{K}}_{\beta }=\left[{\mu }_{\beta }\left(\nabla {\mathbf{v}}_{\beta }+{\left(\nabla {\mathbf{v}}_{\beta }\right)}^T\right)-\frac{2}{3}{\mu }_{\beta }\left(\nabla {\mathbf{v}}_{\beta }\right)\mathbf{I}\right]\mathbf{n}$
Mass transport in gas domain (i − specie = H2, H2S)
$⟨C_{H_2}^{\beta }⟩=C_{H_2}^0;⟨C_{H_{2}S}^{\beta }⟩=C_{H2S}^0=0$ (inlet concentration)	at z = LR	(43)
$-\mathbf{n}\cdot {\mathrm{D}}_i^{\beta }\nabla ⟨C_i^{\beta }⟩=0$ (outlet diffusive contribution)	at z = 0	(44)
$-\mathbf{n}\cdot {\mathbf{N}}_i^{\beta }=0$ (impermeability)	at r = rR and Aωβ	(45)
$-\mathbf{n}\cdot {\mathbf{N}}_i^{\beta }=0$ (symmetry)	at plane in x = 0	(46)
Mass transport in liquid domain (i − specie = R − S, H2, H2S)
$⟨C_{R-S}^{\gamma }⟩=C_{R-S}^0;⟨C_{H_2}^{\gamma }⟩=⟨C_{H_{2}S}^{\gamma }⟩=0$ (inlet concentration)	at z = LR	(47)
$-\mathbf{n}\cdot {\mathrm{D}}_i^{\gamma }\nabla ⟨C_i^{\gamma }⟩=0$ (outlet diffusive contribution)	at z = 0	(48)
$-\mathbf{n}\cdot {\mathbf{N}}_i^{\gamma }=0$ (impermeability)	at r = rR	(49)
$-\mathbf{n}\cdot {\mathbf{N}}_i^{\gamma }=K_{g,i}^{\mathit{\beta \gamma }}\left(⟨C_i^{\beta }⟩R_{g}T/H_i-⟨C_i^{\gamma }⟩\right)$ (exchange flux)	at Aωγ	(50)
$-\mathbf{n}\cdot {\mathbf{N}}_i^{\gamma }=0$ (symmetry)	at plane in x = 0	(51)
Mass transport in solid domain (i − specie = R − S, H2, H2S)
${⟨C_i^{\omega }⟩}^{\omega }=⟨C_i^{\gamma }⟩$ (not mass resistences)	at Aωγ	(52)
$-\mathbf{n}\cdot {\mathbf{N}}_i^{\omega }=0$ (concentration field continuity)	at catalyst centers	(53)
$-\mathbf{n}\cdot {\mathbf{N}}_i^{\omega }=0$ (symmetry)	at plane in x = 0	(54)

Table 4.

Boundary conditions for TBR model.

${\mathbf{N}}_i^j$ is the total flux, given by ${\mathbf{N}}_i^j=-{\mathrm{D}}_i^j\nabla ⟨C_i^j⟩+{\mathbf{v}}_j⟨C_i^j⟩$.

This rigorous model allows to access further than the local fields scale of concentration, velocity, and pressure, such as axial and radial reaction rates, interstitial fluxes for each individual specie, localization of channeling in the flows, and with further post-processing of the data, effectiveness factors for individual pellets, wall-effects analysis, and transfer resistances analysis, to name some examples.

3.5. Kinetic model

In order to give insights into the communication of the information between length scales, a hydrodesulphurization (HDS) reaction for a light gas-oil (LGO) was implemented in both the heterogeneous micropores model and the TBR model. The reaction follows the stoichiometric expression.

$R‐S_{\left(\mathit{liq}\right)}+2H_{2\left(\mathit{gas}\right)}\to R‐H_{2\left(\mathit{liq}\right)}+H_2{S}_{\left(\mathit{gas}\right)}$E55

In this expression, R-S is the sulfurized specie and R-H₂ is the desulfurized specie.

Reaction follows a kinetic expression of type of Langmuir-Hinshelwood/Hougen-Watson [12], described by the following expression.

$-r_{\mathit{HDS}}=\frac{k_{\mathit{HDS}}{⟨C_{R-S}^{\theta }⟩}^{\theta }{\left({⟨C_{H_2}^{\theta }⟩}^{\theta }\right)}^{0.5}}{{\left(1+K_{H_{2}S}{⟨C_{H_{2}S}^{\theta }⟩}^{\theta }\right)}^2};i-\mathit{\text{phase}}=\omega ,\lambda$E56

In this expression, i = ω for the TBR model that considers the catalyst as a pseudo-phase, and i = λ for the heterogeneous micropores model. It should be remembered that it is also required to make the adjustment specified in Eq. (11).

3.6. CFD computation

To solve the CFD models at the three different scales previously described, COMSOL Multiphysics software simulations were implemented. Due to the multiphysic and multiscale nature, as well as to the highly non-linearity and the geometrical complexity of the models; great computational resources, as RAM memory, processing capacity and computational times were required [1]. Thus, the CFD models and representation to study these cases are, in a great extent, limited by computational resources, and therefore important simplification in geometries as symmetry assumptions have been implemented. In spite of the simplifications, great computational resources were necessary, as it can be seen in Table 5.

Model	Mesh$\left[\begin{array}{l}\mathit{\text{Triangular}}\\ \mathit{\text{Elements}}\end{array}\right]$	RAM memory [GB]	Virtual memory [GB]	Computing time [h]
Micropores	6.15 × 10^6	13.6	41.1	0.4
Mass closure vector	2.5 × 10^6	13.4	33.13	0.3
Energy closure vector	6.15 × 10^6	90.84	110.9	0.82
Pseudo-homogeneous	5.6 × 10^4	1.26	1.48	0.004
Reactor model	6 × 10^6	80	110	3–1 week*

Table 5.

Details of the computing resources required in each CFD model.

^*

Depending on the particular simulation/case tested.

The models were solved in a workstation with a dual socket Intel^® Xeon^® E5-2603 v3 processor (15 M of Cache and 1.60 GHz) and 160 GB of RAM memory. It is important to note that a segregated solving method was used in computation with the purpose of obtaining sufficiency in installed RAM memory.

Also, it is important to note that several CFD commercial software have preloaded templates with the commonly used physics in most science and engineering areas. However, for the study and solution of the boundary values problems (BVPs) for the closure vectors there is no existing template, as it is a specific problem of mathematical nature. Then, it was proceeded to approach these problems aided by the “Coefficients Form Partial Differential Equation (PDE)” from COMSOL Multiphysics. The user definitions and variables were also used in order to incorporate the interfacial momentum exchange models that consider the interactions between the three phases, for the reactor model.

3.7. Model assumptions

For the reactor scale model, the TBR is considered to operate in a Trickle regime where the gas and liquid flow descendant at co-current, and that the operation is isothermal. On another hand, the density and viscosity of gas and liquid phases are considered constant, that the catalyst activity do not change with time, and that vaporization and condensation of gasoil do not take place.

In addition, it is assumed that chemical reactions take place only in the solid catalyst, which is considered to be completely wet, for purposes of mass transport model. As for the selection of the value of N_C = D_R/D_P, in order to neglect the wall effects, a value of N_C ≈ 10 is considered sufficient for an accurate prediction of both, pressure drop and holdup. A brief discussion about the results that support the N_C selection is presented within the results section. Also, it is considered that the ordered catalyst bed representation built for the model contains enough characteristics of a real catalytic bed, and that the symmetry assumption can be implemented. Further discussion in the results section support the validity of these assumptions.

For pellet scale model, it is considered that there is no resistance to heat transfer and mass on the surface of the catalyst.

4. Results

4.1. Mass and energy transport at micropores model scale

The solution of boundary value problem of mass and energy transport in porous scale allows the evaluation of concentration fields for R-S, H₂, and H₂S species in the liquid phase that fills the interstitial domain of the catalyst; as well as the temperature fields for the solid matrix and the interstitial fluid. Figure 5 shows the details of the concentrations and temperature fields obtained from the micropores model, described by Eqs. (3)–(10). It is important to point out that the shown fields images were chosen due to the clarity of the scales.

Figure 5.

Concentration fields for the micropores model (a) temperature increase, (b) H₂S generation, (c) H₂ consumption, and (d) dimensionless sulfurated specie concentration.

It is noteworthy that the concentration fields are shown in the entirety of the fluid domain, which is consequence of the consideration that the interstitial fluid domain is completely filled due to capillarity forces. In addition, the model considers that the catalyst is completely wet, implying that the gas phase in the reactor do not contact the solid. Thus, there is no need to model the mass transport with the gas phase. Furthermore, due to these considerations, the gas phase, at this scale, cannot be modelled.

Also, the kinetic expression in Eq. (56) requires the H₂ concentration at the liquid phase; which must be determined by the equilibria between the H₂ concentrations in liquid and gas phases, as well by the mass transport resistances for the hydrogen between liquid and solid phase, Eq. (50). Thus, considering the prior, the H₂, H₂S and the temperature at the catalyst surface is set to be equal to those in the gas phase, Eqs. (6) and (10). These are clearly further simplifications in the modelling of the micropores model. These same simplifications are assumed in the pseudo-homogeneous pellet model.

Regarding the temperature increase field, it is noteworthy that the maximum temperature increase is lower than 4 × 10⁻³ K for a 0.35 mm diameter catalyst, which supports the assumption of an isothermal operation in the reactor model. Also, in the concentration fields, it can be seen that there is a generation of barely more than 1 mol/m³ of H₂S, a consumption of 1.8 mol/m³ of H₂, and a decrease of 6% of the sulfurized specie.

4.2. Mass and heat closure vectors, and effective transport coefficients

Figure 6 shows the closure vectors fields for the mass transfer and heat transfer problems, given by the solution of the boundary values problems of Eqs. (14)–(20). These results allow the evaluation of the effective transport coefficients (D_eff, K_eff), through Eqs. (12) and (13). Those evaluated coefficients are then used in the pseudo-homogeneous catalyst model and the reactor model. These are the bridges that up-scale the geometrical information and characteristics of the porous media in the catalyst microstructure to the upper scales. The tensor components evaluated values are $\left({\epsilon }_{\lambda }{\mathrm{D}}_{\mathit{eff},\mathit{xx}}/{\mathrm{D}}_i^{\lambda }={\epsilon }_{\lambda }{\mathrm{D}}_{\mathit{eff},\mathit{yy}}/{\mathrm{D}}_i^{\lambda }=0.222848\right)$ and (K_eff, xx/k_λ = K_eff, yy/k_λ = 3.84824), it can be seen that both tensors are symmetric.

Figure 6.

Closure vector fields (a) x and y components of the mass transfer closure vector and (b) x and y components of the heat transfer closure vector.

4.3. Mass and energy transport at pseudo-homogeneous catalyst model scale

Figure 7 shows the solution for the boundary values problem specified by Eqs. (21)–(30); for comparison purposes, the selected shown fields are the same than those selected in Figure 5. In order to obtain these fields in the pseudo-homogeneous catalyst model, the evaluated effective transport coefficients shown in the last section were used in the coefficients of Eqs. (21) and (22).

Figure 7.

Concentration fields for the catalyst scale model (a) temperature increase, (b) H₂S generation, (c) H₂ consumption, and (d) dimensionless sulfurated specie concentration.

As it can be seen, both models depict similar behaviors, although in the pseudo-homogeneous catalyst model the HDS reaction seems to be slower compared with the reaction in the micropores model. This difference conduces to an increase of temperature 28% lower, 10% less production of H₂S, 22% less consumption of the H₂, and 1% less decrease in the sulfurized specie concentration. The difference between both models can be attributed to that the geometric information captured by the effective coefficients is not enough, meaning that some information of the micropores structure is not scaling up.

4.4. Reactor model

4.4.1. Validation of the reactor model

It should be clear that, even though the models are theoretically suitable to describe and predict the transport phenomena in the systems, these needs to be validated in order to use their results in further applications with certainty. However, the validation is in many cases quite difficult as there are not enough available experimental data to compare or the experimental information is difficult to acquire, especially when microscales porous media models are intended to be validated, as the punctual microscales phenomena is hardly measurable or observable by experimental methods. Then, models with multiscale modelling approaches such as the presented in this chapter cannot be easily validated in the microscales.

For the micropores scale, the effective transport coefficients evaluated are in the range of the commonly found values in the literature, and thus the geometrical representation of porous structure can be considered as suitable. For the pseudo-homogeneous catalyst model, the mass and energy transport models have been extensively proved and validated in literature, then the model do not require more discussion about its validation.

In the case of the reactor model, the validation was carried out for the hydrodynamic and kinetic behavior. The hydrodynamics were validated against pressure drop and liquid holdup data found in literature [25], achieving mean absolute relative errors (MARE = (∑|Experimental − Predicted|/Experimental)/n), below 5% in the pressure drop prediction and below 8% for liquid holdup. To validate the hydrodynamics, the geometrical characteristics of the bed, as the column and catalyst dimensions and bed porosity, were adapted to be similar to the experimental setup of Al-Dahhan and co-worker [25]. For the kinetic behavior, the reactor geometrical model was adapted to have similar characteristics than reactor described in the theoretical work of Chacón and co-workers [12]. The model exhibits 5.12% difference in the prediction of the sulfurized conversion at 7 mm length (length of the reactor model discussed in this chapter).

4.4.2. Mass transport in the reactor model

In the reactor model scale, the velocities and pressure profiles in the fluid phases are obtained; as well as the concentration profiles for all species in the phases in which they are present. The sulfurized species do not undergo evaporation, and thus it can only be found in the liquid and solid phase; the gas phase is mainly constituted by the hydrogen, and contains a small amount of H₂S that is produced in the solid phase, and is able to transport through the fluid phase and solubilize in the gas phase; Also, the high pressures and temperatures enable the hydrogen to transport to the fluid phase and then to the solid, where the reaction takes place.

Figure 8 depicts concentration fields for the different species at the catalyst domain in the reactor model in a zx cut-plane, in a selected catalyst in the third layer of catalyst from the inlet: (a) dimensionless concentration of the sulfurized specie, (b) H₂S concentration produced by the reaction, (c) H₂ concentration, and (d) liquid to gas H₂ concentrations ratio in the fluid phases in several xy cut-planes in the reactor.

Figure 8.

zx cut-plane in a selected catalyst in the reactor model (a) dimensionless sulfurized specie concentration, (b) H₂S concentration, (c) H₂ concentration, and (d) liquid/gas H₂ concentration ratio in the fluid phases.

It is important to highlight that there are important resistances to the mass transfer of the hydrogen from the gas phase to the liquid phase; so that at most, only 25% of the hydrogen fed in the gas phase transports to the liquid phase. Thus, the hydrogen that reaches the catalysts is considerably less than the quantity considered for the micropores and pseudo-homogeneous models. In Figure 8(d), it can also be seen that there is a wall-effect of around two catalyst diameters.

4.5. Comparison between scales

Figure 9 shows the H₂ concentration fields obtained in the different scales models: (a) micropores model, (b) pseudo-homogeneous model, (c) selected catalyst, in the third catalyst layer from the inlet, in the reactor model, and (d) shows the liquid to gas H₂ concentrations ratio, as well as the selected catalyst.

Figure 9.

H₂ concentration fields at (a) micropores scale, (b) pseudo-homogeneous catalyst scale, (c) selected catalyst at reactor scale, and (d) selected catalyst particle with liquid/gas H₂ concentration ratio.

As it can be seen, in the micropores and the pseudo-homogeneous catalysts models, the H₂ concentration fields inside the catalyst particle domain are very similar, and show a symmetric distribution of the field, as previously discussed. However, this symmetrical behavior is not observed in the catalyst particle in the reactor model, and the H₂ concentration variations are in the range of 59.5–121.9 mol/m³, which is a greater range that the one observed in the lower scales. And however, as seen by Figures 5, 7 and 8, there is a lower consumption of the sulfurized specie in the selected catalyst in the reactor model. This is due to the interactions between the mass transport inside the catalysts, and the mass transport and hydrodynamics inside the reactor. Considering the prior observations, these exhibit that there is a need to stablish multiscale analyses and bridges that communicate the up-scales phenomena to the down-scales in order to properly study the mass transport phenomena occurring in the micropores scales.

Another important aspect to highlight in Figure 9(c) is that a lower concentration can be seen at the top and higher concentration at the bottom of the catalyst. It is particularly noteworthy since the H₂ concentrations are expected to decrease when moving down in the reactor due to the reaction. The reason for this behavior can be found in Figure 9(d), from which it can be deduced that much of the distribution/transport of the H₂ in the catalysts due to the distribution and transport of this specie in the liquid phase; which can be seen that increases when moving down in the reactor. In addition, it should be reminded that it was assumed that there are no resistances in the transport of the hydrogen from the liquid to the solid phase, which is specified by Eq. (56).

Further information on the H₂ concentration field behavior can be seen in Figure 10, where the total H₂ fluxes are shown inside the selected catalyst and in the surrounding interstitial fluid. These flux arrows are pointing the same direction in which the hydrogen is being dragged. In addition, it is observed that a great part of the H₂ behavior inside the catalyst is due to the concentration gradients inside the catalyst, which originate diffusive transport, and then it is not only due to the consumption for the HDS reaction. It can also be seen that there is an important influence of the surrounding fluid, which has diffusive and convective contributions to the mass transport, being the later greater than the diffusive.

Figure 10.

H₂ fluxes inside the selected catalyst and its surroundings fluid.

It is important to note that this kind of analyses are possible as a result of the comparison of the information obtained from the solution of the different length scales models (multiscale analysis), as well as the establishment of detailed mathematical models for the phenomena involved (multiphysics analysis).

5. Conclusions

Simulation and analyses of an HDS process at different scales were carried out; the micropores and pseudo-homogeneous catalyst models included three mass transport equations for the consumed and generated species, two heat transfer equations for the micropores model, an one heat transfer equation for the pseudo-homogeneous catalysts model; at the reactor-scale, the model consisted of two momentum balance equations, and eight mass transport equations, for the consumed and produced species and their transport between phases. These models allow the multiphysics and multiscale analyses of the process, which is possible as the different length scales share information between them.

In the micropores representation model, with a 2:1 scale, the study and determination of the effective transport coefficients for the mass and heat transfer was carried out. The geometrical representation was developed from a vectorized micrograph of a real porous media obtained from literature. The evaluated coefficients are in the order of the ones found in literature, and can be thus considered as suitable.

In literature, it considered that the solution of the mass and heat transfer with superficial reaction at a micropores scale is hardly developed, and then average mass and heat transport equations are used (pseudo-homogeneous models). The proper equivalence between these two approaches is achieved through adequate values for effective transport coefficients. In the analyses shown here, it can be seen that similar tendencies in the concentration and temperature fields are achieved, but also differences that may be significant are shown.

The CFD reactor model for the HDS process has a good agreement with the experimental pressure drop and liquid holdup data, and with the theoretical conversion data. Thus, the model is considered as suitable for the prediction of hydrodynamics and kinetics behaviors, and can be used for further analyses.

The H₂ mass transport with reaction analysis in a catalytic particle at the different scales models allow to show that the upper scales phenomena, as the effect of the mass transport in the interstitial bed fluid, modifies in a great extent the behavior of the mass transport in the catalytic particle; that the micropores and the pseudo-homogeneous models are unable to represent that behavior; and that only a multiscale analysis allows to study and analyze these kinds of phenomena.

Also, regarding the implementation of the models in the CFD commercial software COMSOL Multiphysics, it was necessary to specify the boundary values problem (BVP) for the closure vector through the general COMSOL interphase for the “Coefficients Form PDE”, suitable for many well-known PDEs systems; and for the reactor model, it was necessary to use the user definitions and variables to specify the interfacial momentum exchange models.

The three models at the different scales have been specified by commonly found boundary conditions for BVP with closed domains, this is Dirichlet, Neumann and Robinson boundary conditions. In order to specify that some information of the field, the field derivate or a combination of the field and its derivate, is known.

With the CFD reactor model, it is possible to analyze the wall effect, convective and diffusive fluxes for all species in every phase, reaction rate, and with further postprocessing, effectiveness factors, which allows to determine the effect of the T, P operation conditions, as well as the gas and liquid inflow velocities over the reactor conversion.

Nomenclatures

Abbreviations
MP	micropores
PHC	pseudo-homogeneous catalyst
RS	reactor scale
HDS	hydrodesulphurization
TBR	trickle bed reactor
dp	micropore diameter
Dp	Pellet diameter
rR	reactor radius
LR	reactor length
Aij	interfacial area between phase i and j
Dij	diffusion of specie i in the phase j
kj	specie j conductivity
Nc	critical reactor to pellet diameters ratio
0	at initial conditions
σ	solid
β	gas
γ	liquid
λ	catalyst micropores fluid
Λ	catalyst micropores fluid
ω	pseudo-homogeneous solid
εi	volume fraction of phase i