Computational Flow Modelling of Multiphase Reacting Flow in Trickle-bed Reactors with Applications to the Catalytic Abatement of Liquid Pollutants

1. Introduction

The contemporary status of computational flow modelling encourages strongly the application of modern CFD codes on the design and investigation of process equipment and aims to accomplish novel and stringent environmental regulations on the decontamination of high strength wastewaters. The application of CFD is envisaged as a prevailing tool to gain more control in the delivery of reactants or energy and in the removal of products since any chemical or physical transformation process requires the addition or removal of different materials and energy. Therefore, the distribution of materials as well as the energy within the process vessel is improved noticeably if one has the capability to predict and control fluid dynamics by means of CFD. Instead of extensive and costly experimentation, CFD models allow the optimization of process equipment to maximize industrial benefit following rigorous environmental regulations. CFD enable the establishment of quantitative relationships between the throughput and operating flow regime avoiding a massive number of empirical parameters. Transport parameters such as heat- mass-transfer coefficients and several hydrodynamic parameters including film structures and thicknesses, boundary layer geometries, pressure drop, wetted area or even the residence time distribution are investigated and optimized for multiple process configurations. Consequently, CFD aids in the identification of most appropriate process configurations from the large selection of alternatives.

In fact, in the realm of environmental reaction engineering, major improvements on active and stable heterogeneous catalyst are intended to be realized in the near future. This fact will allow substantial enhancements of catalytic wet oxidation technology to become economically and environmentally feasible in comparison to most of the conventional treatments. In this ambit, agro-food processing wastewaters will be undertaken as an example of non-biodegradable and phytotoxic effluents. Polyphenolic compounds deserve particular attention in what concerns the process performance at different levels as follows: total organic carbon reduction, intermediate compounds formation, catalyst stability and kinetic modelling. The kinetic parameters will be obtained taking into account the generalized lumped kinetic model and this information will be further integrated into a hydrodynamic model developed by means of CFD where the reaction aspects as well as the transport mechanisms were accounted for a high-pressure trickle-bed reactor. Two multiphase CFD frameworks – Euler-Euler and Volume-of-Fluid – will be evaluated for either cold flow or reacting flow conditions. Axial and radial profiles of hydrodynamic and reaction parameters will be thoroughly evaluated at different operating conditions. The influence of gas and liquid flow rates as well as the effect of temperature and pressure will be investigated in terms of total organic carbon conversions.

2. Literature review

In the field of environmental reaction engineering, long-established technologies have been intensified during the last decade to demonstrate new potentialities and hindrances from a process viewpoint. Having recognized the effect of technological threats generated by dangerous substances on the environment, it is possible to identify a gap of knowledge on how to provide consistent remediation technologies of contaminated air, soil and water. In this regard, multiphase reactors such as trickle-bed reactors (TBRs) can then be evaluated on its further deployment for the catalytic abatement of moderate- to high-strength wastewaters if one takes into account modern and advanced simulation tools as CFD codes.

The specialized literature identified various dynamic circumstances in TBRs as a new breed potential for performance enhancement arising from the interaction between the gas and liquid phases in supplying reactants to the catalyst surface and pores. This fact is even reinforced due to the multifaceted relationship between fluid dynamics, catalyst wetting and chemical reaction. Multiphase reactors such as TBRs have been commonly designed and operated under steady-state mode in the industrial applications. However, novel operating modes under transient conditions has generated an emergent interest in view of the significant performance enhancement encompassing chemical conversion and products selectivity (Dudukovic et al., 1999; Boelhouwer et al., 2002; Silveston & Hanika, 2002, 2004). The transient operation of TBR revealed key characteristics reinforced by the dynamic nature of gas-liquid flow over the catalyst packing, which ultimately leads to the cumbersome non-linearity behaviour of trickle beds (Khadilkar et al., 2005). Here, we pursue our prior efforts to provide new insights on environmentally-based TBRs to detoxify hazardous polluted liquid streams. This task is sought to be accomplished by the experimental work devoted to investigate the suitability of this methodology as well as by state-of-the-art simulation tools considered by CFD codes. The case study encompasses the application of a trickle-bed reactor to the detoxification of phenolic wastewaters through a series of experimental runs accompanied and validated by a multiphase Volume-of-Fluid (VOF) model formerly reported (Lopes & Quinta-Ferreira, 2009). At a first glance, a brief prospect of literature on trickle-bed reactor modelling and validation is presented at reaction conditions. Then, the application of CFD is reinforced by its investigation on the minimization of hazardous materials.

Iliuta et al. (2006) investigated how plugging affect hydrotreating performance by simulating simultaneous fines deposition under catalytic in hydrodesulfurization in TBRs. The three-phase heterogeneous model developed to simulate the trickle bed performance incorporates the concentration gradients inside the catalyst particle and solid deposit. The model correctly handled the effect of the fine particles deposition process being identified the bed plugging and the increase of the resistance to two-phase flow. Ayude et al. (2005) presented a phenomenological approach to interpret the effect of liquid flow modulation in TBRs at the particle scale. Dynamic reactant profiles inside the catalytic particle were obtained for different cycling and system conditions being the experimental trends explained with this approach. Liu et al. (2005) examined the periodically operated TBR for ethylanthraquinones hydrogenation. A dynamic model was verified by simulating the effect of cycle period and split on the conversion and the selectivity enhancement and compared with the experimental results. An experimental and theoretical study of forced unsteady-state operation of TBR was presented by Lange et al. (2004) in comparison to the steady-state operation. The simulation studies demonstrate that the liquid flow variation has a strong influence on the liquid hold-up oscillation and on the catalyst wetting efficiency. Rajashekharam et al. (1998) developed a TBR model for the hydrogenation of 2,4 dinitrotoluene with experimental verification. The model incorporates the contributions of partial wetting and stagnant liquid hold-up effects in addition to external and intraparticle mass transfer resistances for a complex consecutive/parallel reaction scheme under consideration represented by Langmuir-Hinshelwood kinetics.

On the investigation of CFD for the risk assessments related to industrial environments, Calay & Holdo (2008) have modelled the dispersion of flashing jets in where gas is kept in liquid form under high-pressure. The authors presented reproducible results to those phenomena related to turbulence, droplet transport, evaporation, break-up and coalescence events. Gavelli et al. (2010) presented a reliable method to quantify and bound LNG vapour dispersion hazard distances for regulatory purposes. The mathematical models encompass conservative assumptions in the derivation of the LNG vapour source term. It is suggested that sophisticated models could be used to predict the flow and vaporization of LNG by numerical integration of the shallow water equations. Pontiggia et al. (2010) have assessed the consequences of hazardous gas releases in urban areas through CFD simulations. Once again, CFD has demonstrated its versatility to represent realistically such complex geometries in detailed simulation models leading to three-dimensional flow fields that strongly influence gas dispersion. The authors have pointed out that integral models can both overestimate and underestimate the magnitude of consequences related to hazardous material releases in urban areas.

From the above survey, CFD models can be trustfully applied in different frameworks, nevertheless such approach has scarcely been reported to predict the reaction behaviour of trickle beds. Apart from the hydrodynamic models evolved at cold flow conditions, we extend our previous VOF model to gain further insight on how process parameters such as liquid velocity, surface tension and wetting phenomena affect the multiphase flow in high-pressure trickle-bed reactor. A Volume-of-Fluid model was developed to simulate the multiphase reactive flow in the catalytic wet oxidation of mimicked agro-industry effluents. As long as the temperature and liquid flow rate play a major role on the decontamination rate of liquid pollutants, several computational runs have to be validated accordingly and performed under transient conditions. This methodology provides a more fundamental understanding of flow environment and TBR chemical behaviour.

3. Mathematical models

3.1. Governing flow equations

A trickle bed based on a cylindrical geometry (50 mm_ID×1.0 m_Length) was modelled with a specified void fraction and a set of fluid physical properties. The computational geometry was designed so that a distance gap of about 3% of the sphere diameter facilitates the grid generation avoiding numerical difficulties that arise in the calculation of convective terms as described elsewhere (Lopes & Quinta-Ferreira, 2008; Nijemeisland & Dixon, 2001). The VOF method was used to compute velocity field, liquid volume fraction distributions as well as the total organic carbon concentration that was used to quantify the mineralization degree of liquid pollutants. The multiphase flow is assumed to be vertical downward and incompressible, with the mathematical description based on the Navier-Stokes equations for momentum and mass conservation of a viscous fluid flow through a three-dimensional catalytic bed The variable fields are shared by both phases and correspond to volume-averaged values knowing the volume fraction, α_q, of each phase, q, in the entire computational domain. The volume fraction equation for the phase q is given by:

$$\frac{\partial }{\partial t}\left({\alpha }_q{\rho }_q\right)+\nabla \cdot \left({\alpha }_q{\rho }_q{U}_q\right)=0\text{with}q=g\text{or}l$$E1

where g and l denote, respectively, the gas and liquid phases, t, being the time, and through the resolution of the momentum equation shared by the two considered fluids:

$$\frac{\partial }{\partial t}\left({\alpha }_q{\rho }_q{U}_q\right)+\nabla \cdot \left({\alpha }_q{\rho }_q{U}_q{U}_q\right)=-{\alpha }_q\nabla p+{\alpha }_q{\rho }_{q}g+\nabla \cdot {\alpha }_q\left(\overline{\overline{{\tau }_q}}+\overline{\overline{{\tau }_{t,q}}}\right)+I_q\text{with}q=g\text{or}l$$E2

where p, g and the physical properties (density, ρ; and viscosity, μ) being determined by volume-weighted averages. I_q is the interphase momentum exchange term and ${\overline{\overline{\tau }}}_q$ and ${\overline{\overline{\tau }}}_{t,q}$are, respectively, the viscous stress tensor and the turbulent stress tensor, defined as follows:

$$\overline{\overline{{\tau }_q}}={\mu }_q\left(\nabla U_q+\nabla U_q^t\right)+\left({\lambda }_q-\frac{2}{3}{\mu }_q\right)\nabla U_q\overline{\overline{I}}$$E3

and

$${\overline{\overline{\tau }}}_{t,q}={\mu }_{t,q}\left(\nabla U_q+\nabla U_q^t\right)-\frac{2}{3}\left(k_q+{\mu }_{t,q}\nabla U_q\right)\overline{\overline{I}}$$E4

3.2. Free surface model: surface tension and wall adhesion

The continuum surface force model proposed by Brackbill et al. (1992) has been used to compute the surface tension. The pressure drop across the surface depends upon the surface tension coefficient, σ, and the surface curvature as measured by two radii in orthogonal directions, R₁ and R₂, as expressed by Eq. 5.

$$p_2-p_1=\sigma \left(\frac{1}{R_1}+\frac{1}{R_2}\right)$$E5

where p₁ and p₂ are the pressures in the two fluids on either side of the interface. The surface curvature is computed from local gradients in the surface normal at the interface. n is the surface normal, defined as the gradient of α_i :$n=\nabla {\alpha }_i$. The curvature, κ, is defined in terms of the divergence of the unit normal,$$\widehat{n}$$: $$\kappa =\nabla \cdot \widehat{n}$$where$$\widehat{n}=n/\left|n\right|$$. The forces at the surface are expressed as a volume force using the divergence theorem assuming the form of Eq. 6.

$$F_j={\displaystyle \sum _{pairsij,i<j}{\sigma }_{ij}\frac{{\alpha }_i{\rho }_i{\kappa }_j\nabla {\alpha }_j+{\alpha }_j{\rho }_j{\kappa }_i\nabla {\alpha }_i}{\frac{1}{2}\left({\rho }_i+{\rho }_j\right)}}$$E6

3.3. Species continuity and energy equations

The predicted flow field including velocities and volume fractions of both phases was further used to solve species transport equations in the catalytic wet air oxidation of a model phenolic solution in the trickle-bed reactor. These equations are expressed in the mass balance equation for any species, i:

$$\frac{\partial {\alpha }_q{\rho }_q{C}_{q,i}}{\partial t}+\nabla \left({\alpha }_q{\rho }_q{u}_q{C}_{q,i}\right)=\nabla \left({\alpha }_q{\rho }_q{D}_{i,m}\nabla C_{q,i}\right)+{\alpha }_q{\rho }_q{S}_{q,i}$$E7

where, C_q,i is the concentration of species i in the q^th phase (gas or liquid), ρ_q and α_q is the density and volume fraction of the q^th phase. S_q,i is the source for species i in phase q. Volume averaged properties of fluids were used for calculating the flux across the control cell. Two-film theory was used for accounting mass transfer and the resistance in gas-liquid film was considered as the rate limiting resistance (Bhaskar et al., 2004). Mass transfer coefficient was computed according to the Satterfield et al. (1978) correlation and heat transfer coefficient was calculated according to the correlation developed by Boelhouwer et al. (2001) as expressed by equations 8 and 9, respectively.

$$Sh=0.\text{815}R{e}^{0.\text{822}}S{c}^{\text{1}/\text{3}}$$E8

$$Nu=0.\text{111}R{e}^{0.\text{8}}P{r}^{\text{1}/\text{3}}$$E9

The energy equation, also shared among the phases, is shown in equation 10:

$$\frac{\partial }{\partial t}\left(\rho E\right)+\nabla \cdot \left(\overrightarrow{v}\left(\rho E+p\right)\right)=\nabla \cdot \left(k_{eff}\nabla T\right)+S_h$$E10

The VOF model treats energy, E, and temperature, T, as mass-averaged variables:

$$E=\frac{{\displaystyle \sum _{q=1}^n{\alpha }_q{\rho }_q{E}_q}}{{\displaystyle \sum _{q=1}^n{\alpha }_q{\rho }_q}}$$E11

where E_q for each phase is based on the specific heat of that phase and the shared temperature. The properties ρ and k_eff (effective thermal conductivity) are shared by the phases. The source term, S_h, contains contributions from volumetric reaction heat sources, which is given by the product ΔH×ρ×r_TOC.

3.4. Two-phase k-ε turbulence model

Taking into account that the Reynolds numbers range for the gas phase is wide (min: 10, max: 2500), the mixture k-ε approach is used for turbulence modelling (Elghobashi et al., 1984). For incompressible flows, the turbulence parameters are calculated from equations 10-11:

$$\frac{\partial }{\partial t}\left({\rho }_m\epsilon \right)+\nabla \cdot \left({\rho }_m{u^{\prime }}_{m}k\right)=\nabla \cdot \left(\frac{{\mu }_{t,m}}{{\sigma }_k}\nabla k\right)+G_{k,m}-{\rho }_m\epsilon$$E12

$$\frac{\partial }{\partial t}\left({\rho }_m\epsilon \right)+\nabla \cdot \left({\rho }_m{\overrightarrow{u}}_m\epsilon \right)=\nabla \cdot \left(\frac{{\mu }_{t,m}}{{\sigma }_k}\nabla \epsilon \right)+\frac{{\epsilon }_m}{k_m}\times \left(C_{1\epsilon }{G}_{k,m}-C_{2\epsilon }{\rho }_m\epsilon \right)$$E13

and the turbulent viscosity μ_t,m and the production of turbulence kinetic energy, G_k,m are computed from Eq. 14.

$$\begin{array}{l}{\mu }_{t,q}={\rho }_q{C}_{\mu }\frac{k_q^2}{{\epsilon }_q}\\ G_{k,m}={\mu }_{t,m}\left(\nabla {\overrightarrow{u}}_m+{\left(\nabla {\overrightarrow{u}}_m\right)}^T\right):\nabla {\overrightarrow{u}}_m\end{array}$$E14

C_µ is equal to 0.09 and C_1ε and C_2ε are the constants of standard k-ε model: 1.44 and 1.92, respectively, whereas σ_k and σ_ε are the turbulent Prandtl numbers for k and ε, 1.0 and 1.3, respectively.

3.5. Numerical simulation

The interstitial space of the trickle-bed reactor was designed through a tetrahedral mesh representing that was created using the integrated solid modelling and meshing program GAMBIT (2005) mimicking the characteristic dimensions of commercial catalyst N-140 supplied by the Süd-Chemie Group, Munich. The VOF method simulates free-surface flow by means of a fluid fraction function, which has a value between unity and zero. The discretization of the governing equations is done by the finite-volume method. The grid independency was established after the evaluation of different mesh natures and apertures in order to isolate mesh related discretization errors. All transport equations were discretized to be at least second order accurate in space. A segregated implicit solver available in commercial CFD package FLUENT 6 (2005) was employed to evaluate the resulting linear system of equations. The conditions required for grid convergent results are based on a 1% relative error criterion and the simulations accuracy has been assessed by comparisons to experimental data available in the literature. At the interface, the additional interaction conditions depend on interfacial velocity and gradient of the surface tension.

The CWAO kinetic parameters for the commercial catalyst N-140 were similarly derived to the work developed by Lopes et al. (2007). The right-hand side term of Eq. 7, S_i,q, include the reaction rates in terms of the total organic carbon concentration of the lumped species A, B and C as represented by Eq. 15:

$$\begin{array}{l}-r_{TO{C}_A}=-\frac{d{C}_{TO{C}_A}}{dt}=\left({k^{\prime }}_1+{k^{\prime }}_2\right)C_{TO{C}_A}\\ -r_{TO{C}_B}=-\frac{d{C}_{TO{C}_B}}{dt}={k^{\prime }}_3{C}_{TO{C}_B}-{k^{\prime }}_2{C}_{TO{C}_A}\end{array}$$E15

where first order reactions were assumed for each mechanism step of the Generalized Kinetic Model. GKM considers three types of compounds: easier degraded reactants (A); intermediates with difficult degradation (B) and desired end products, namely carbon dioxide and water (C). In the oxidation process of the phenolic solutions with N-140 catalyst, phenol and acetic acid were formed as intermediate compounds, being totally degraded during the treatment and the overall TOC practically reduced to zero. After integrating these equations a mathematical expression for TOC evolution is obtained in Eq 16:

$$\frac{C_{TOC}}{C_{TO{C}_0}}=\frac{{k^{\prime }}_2}{{k^{\prime }}_1+{k^{\prime }}_2-{k^{\prime }}_3}{e}^{-{k^{\prime }}_{3}t}+\frac{{k^{\prime }}_1-{k^{\prime }}_3}{{k^{\prime }}_1+{k^{\prime }}_2-{k^{\prime }}_3}{e}^{-\left({k^{\prime }}_1-{k^{\prime }}_2\right)t}$$E16

The activation energies and the pre-exponential factors were calculated by using the Arrhenius plot for the N-140 kinetic studies. These values were used in the corresponding expressions of the reaction rate constants k^’₁, k^’₂^, k^’₃ as a function of temperature, according to Arrhenius law as described in Eq 17:

$$\begin{array}{l}{k^{\prime }}_1=452\cdot \exp \left(-\frac{3.121\times {10}^3}{T}\right){\min }^{-1};\\ {k^{\prime }}_2=28.1\cdot \exp \left(-\frac{3.612\times {10}^3}{T}\right){\min }^{-1};\\ {k^{\prime }}_3=4.32\times {10}^6\cdot \exp \left(-\frac{9.814\times {10}^3}{T}\right){\min }^{-1}\end{array}$$E17

The trickle-bed reactor studies have been carried out in pilot plant as illustrated in Fig. 1 comprising a cylindrical reactor in stainless steel (SS-316) with 50 mm of internal diameter and 1.0 m length and the experimental procedure has been described elsewhere (Lopes & Quinta-Ferreira, 2010).

Figure 1.

Schematic diagram of trickle-bed reactor experimental set-up

4. Results and discussion

4.1. Mesh and time sensitivity tests

To check grid independence, several computational runs were carried out to figure out the mesh density with different apertures. An introductory parametric optimization was also accomplished querying the effect of different differencing schemes on the momentum and volume fraction balance equations. High-order differencing schemes based on Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM) and High Resolution Interface Capturing (HRIC) schemes were found to agree better with the experimental data from the literature given that its formulation includes inherently the minimization of artificial numerical dissipation as described elsewhere (Lopes & Quinta-Ferreira, 2009).

We begin with a base case querying the influence of numerical solution parameters including different mesh apertures. The refinement is performed in eight levels: four levels for the coarser meshes with 2×10⁵, 6×10⁵, 1×10⁶, 1.4×10⁶ of tetrahedral cells and four levels for the finer meshes with 1.8×10⁶, 2.2×10⁶, 2.6×10⁶ and 3×10⁶ of tetrahedral cells. Figure 2 displays the influence of the number of tetrahedral cells on the total organic carbon concentration profile when the catalytic abatement of phenolic wastewaters was simulated at L=6 kg/m²s, G=0.3 kg/m²s and T=200 ºC, P=30 bar.

Figure 2.

Comparison of CFD predictions on normalized Total Organic Carbon concentration for different mesh resolutions (L = 6 kg/m²s, G = 0.3 kg/m²s, P = 30 bar, d_p=2 mm)

As can be seen from Fig. 2, the increase of mesh density led to an asymptotic solution as one increases the number of cells from 2.2×10⁶ onwards. Low mesh density characterized by 2×10⁵ of tetrahedral cells at particle surface led to erroneous solutions due to an incorrect definition of boundary layer. As long as the mesh density increases, the theoretical predictions of total organic carbon conversion improve noticeably. Concerning the mesh sensitivity analysis, several computational runs were additionally performed changing the mesh density in the catalyst particle surface to properly capture the boundary layer now on the temperature profile. The effect of mesh aperture on the thermal behaviour of the trickle-bed reactor is portrayed in Fig. 3 at L=6 kg/m²s, G=0.3 kg/m²s and T=200 ºC, P=30 bar.

Figure 3.

Comparison of CFD predictions on normalized bulk temperature for different mesh resolutions (L = 6 kg/m²s, G = 0.3 kg/m²s, P = 30 bar, d_p=2 mm)

Once again, the CFD predictions converged as soon as we used a computational grid comprising 2.2 million of tetrahedral cells. Another mathematical output information revealed by Figs. 2 and 3 is related to the monotony of the total organic carbon conversion and temperature profiles. As evidenced by the mesh sensitivity analysis, the asymptotic demeanour of total organic carbon profile is monotonically decreasing until the final conversion approaches the experimental value at steady-state. On the contrary, the temperature profile is a monotonically decreasing function as one increases the density of the computational grid. Provided that the catalytic wet oxidation follows the Arrhenius's law, the higher conversions or lower percentages of the ratio C_TOC/C_TOC0 are accompanied with increasing values of the bulk phase temperature, see Fig. 3.

The numerical integration encompassed by the CFD model is also dominated by the time-stepping methodology. In the CFD multiphase calculations, we employed first-order and second-order time discretization methods. This latter scheme was found to give confident computed results in comparison to the former one as long as the second-order time-accurate scheme achieves its accuracy by using an Euler backward approximation in time. This fact can in all likelihood be ascribed to the fully implicit methodology which highlighted an unconditional stability. A fair implementation of one of these methods for solving the scalar transport equation entails more than the time-stepping formula. Additionally, we have investigated the time step in the range: 10^-5, 10^-4, 10^-3 and 10^-2 s. Figure 4 plots the CFD computations of total organic carbon conversions attained with those values at L=6 kg/m²s, G=0.3 kg/m²s and T=200 ºC, P=30 bar.

Figure 4.

Effect of time step on normalized Total Organic Carbon concentration (L = 6 kg/m²s, G = 0.3 kg/m²s, P = 30 bar, d_p=2 mm)

A similar asymptotic trend was depicted in Fig. 4 so we have identified a time step of 10^-5 s which agreed better with the experimental data. Concomitantly, the analysis of the temperature profile obtained at L=6 kg/m²s, G=0.53kg/m²s and T=200 ºC, P=30 bar reinforced the above-mentioned selection on the time step assessment, see Fig. 5.

Figure 5.

Effect of time step on normalized bulk temperature (L = 6 kg/m²s, G = 0.3 kg/m²s, P = 30 bar, d_p=2 mm)

It is worthy to notice the reverse behaviour that featured the total organic carbon concentration and temperature profiles illustrated in Figs. 4 and 5. As discussed earlier on the monotonic nature of mesh density functions, this suggests again the Arrhenius's temperature dependence of the rate constant for the chemical oxidation reaction.

4.2. Effect of temperature

At L=2 kg/m²s, G=0.3kg/m²s, P=30 bar, the VOF computations of total organic carbon conversion along the axial coordinate of trickle-bed reactor are shown with experimental data at different temperatures in Fig. 6.

Figure 6.

Mean bulk total organic carbon conversion profiles for axial coordinate at different temperatures (L = 2 kg/m²s, G = 0.3 kg/m²s, P = 30 bar, d_p=2 mm)

In the range 160-200ºC, five characteristic temperatures of catalytic wet oxidation have been used to query the effect of this reaction parameter on organic matter decontamination. As expected, conversion is favoured by increasing the temperature, hence we have the complete mineralization of organic matter after one half of reactor length when the simulation was carried out at the highest temperature (200ºC), see Fig. 6e. However, none experimental measure on total organic carbon concentration has indicated that we have an absolute depletion of organic pollutants from the liquid phase. Examination of successive plots of Fig. 6, provides the ability to compare quantitatively CFD model predictions with the corresponding experimental conversion either through the averaged relative or absolute errors. At the lowest temperature (T=160ºC, see Fig. 6a), the averaged relative errors obtained at z=0.2, 0.4, 0.6, 0.8 and 1m were 16.2, 24.6, 43.1, 46.4 and 44.8%, respectively. Despite these values are remarkably increasing, the averaged absolute errors are decreasing and we had achieved 5.9, 3.6, 3.8, 2.9, 2.1 for z=0.2, 0.4, 0.6, 0.8, 1m, respectively. In fact, the magnitude of the relative errors can be explained by the narrow range of TOC concentrations that were obtained at the reactor egress. At z=0.6 and 0.8m the TOC conversions were higher than 90% which produced very low concentration ratios C_TOC/C_TOC0 as illustrated in Fig. 6. This behaviour on the complementary nature of both averaged relative and absolute errors warned that every validation accomplished with VOF multiphase model should follow a systematic comparison between those rating parameters. Indeed, at T=170ºC the averaged relative errors obtained with CFD model were 18.6, 30.6, 57.8, 64.1 and 65.9% for z=0.2, 0.4, 0.6, 0.8 and 1m whereas the averaged absolute errors were 6.4, 3.8, 4.1, 2.9, 2.1, respectively. Moreover, at the highest temperature the magnitude of the averaged relative errors was enlarged hastily to 32.5, 78.4, 90.6, 90.9 and 91.4% albeit the averaged absolute errors kept almost the same intensity that was early identified at lower temperatures as 4.1, 3.9, 3.8, 2.9 and 2.4 for z=0.2, 0.4, 0.6, 0.8 and 1m, respectively.

According to all VOF computations that have been carried out at L=2 kg/m²s, G=0.3kg/m²s, P=30 bar, one was able to identify the recurring overestimation of total organic carbon conversion for every simulated temperature. Generally, it can be seen from Fig. 6 that computed results agreed better with experimental data at the lower temperatures. This fact has been explained by the formulation of enthalpy balance meaning that a separate enthalpy equation for each phase in Eulerian models is more efficient on the prediction of temperature, which is used to compute the species concentration, against the shared formulation for enthalpy and temperature variables accounted by the current VOF model.

At L=2 kg/m²s, G=0.3kg/m²s, P=30 bar, the averaged axial temperature profiles are depicted by Fig. 7. The thermal behaviour exhibited by the trickle-bed reactor is a distinguishing feature of exothermic reactions. As long as the catalytic wet air oxidation of phenolic pollutants releases energy available from the lower enthalpy of reaction products in comparison to the higher enthalpy of reactants, the chemical conversion of those compounds into water and carbon dioxide is responsible by the temperature increase of bulk phase. This fact led to the heating of catalyst particles more pronounced in the first one fourth of trickle-bed length as portrayed in the thermal profiles of Fig. 7. The higher the inlet temperature was the stronger the increase of bulk phase temperature was attained at steady-state. The averaged relative errors at the lowest temperature were -1.7, -0.5, -0.4, -0.4 and -0.5% for z=0.2, 0.4, 0.6, 0.8 and 1m, respectively. These averaged relative errors were almost identical to the ones obtained at T=170ºC and they have become lower when the catalytic wet oxidation was performed at the highest temperature giving -1.1, -0.5, -0.6, -0.5, and -0.3% for z=0.2, 0.4, 0.6, 0.8 and 1m, respectively.

Figure 7.

Mean bulk temperature profiles for axial coordinate at different temperatures (L = 2 kg/m²s, G = 0.3 kg/m²s, P = 30 bar, d_p=2 mm)

Having analysed quantitatively these findings, it can be argued that the VOF computations of total organic carbon conversion agreed better with experimental data at lower temperatures, whereas the computed VOF predictions of temperature profiles managed better at higher ones.

4.3. Effect of liquid flow rate

The VOF multiphase model has also been used to query the effect of liquid flow rate on total organic carbon conversion within the above-mentioned temperature range, 160-200ºC. Figure 8 shows the axial profiles of total organic carbon conversion attained theoretically and experimentally at L=6 kg/m²s, G=0.3kg/m²s, P=30 bar.

Figure 8.

Mean bulk total organic carbon conversion profiles for axial coordinate at different temperatures (L = 6 kg/m²s, G = 0.3 kg/m²s, P = 30 bar, d_p=2 mm)

A fivefold increase on liquid flow rate led to lower decontamination rates as it can be observed from Fig. 8. This behaviour can be explained by the lower residence time of liquid phase which limited the mineralization degree expressed by reminiscent total organic carbon concentration. The maximum conversion was obtained at the highest temperature giving 90.7% while 85.1% was achieved at T=160ºC. Once again, the higher the inlet temperature is, the improved the total organic carbon conversion become at steady-state. When one increases the liquid flow rate up to L=6 kg/m²s, the VOF multiphase model also overestimated the decontamination rate as already categorized at lower liquid flow rates. The five series of simulated temperatures illustrated in Fig. 8 were additionally employed on the calculation of both averaged relative and absolute errors. At the lowest temperature, the averaged relative errors attained at z=0.2, 0.4, 0.6, 0.8 and 1m were 18.2, 25.4, 32.1, 30.9 and 30.2%, whereas the averaged absolute errors were 10.1, 10.7, 10.8, 9.1 and 5.9, respectively, see Fig. 8a. Increasing the temperature up to T=180ºC, the averaged relative errors became 20.7, 27.8, 33.2, 29.6 and 23.7% while the averaged absolute errors turned out to be 10.9, 10.3, 8.7, 5.1 and 3.3 for z=0.2, 0.4, 0.6, 0.8 and 1m, respectively. If one increases even further the reactor temperature up to T=200ºC, the computed averaged relative errors were 24.7, 31.3, 35.7, 29.5 and 18.4 whereas the averaged absolute errors happened to be 11.8, 10.7, 8.8, 4.7 and 2.1, see Fig. 8e. This in return translates into an unbiased trend of deviations between the computed results and experimental data in dissimilarity with previous results obtained at the lower liquid flow rate, L=2 kg/m²s, keeping invariable the gas flow rate and operating pressure (G=0.3kg/m²s, P=30 bar).

Figure 9 shows the computed results and experimental data of the temperature profiles attained axially when the catalytic wet oxidation has been performed at the same five-valued temperature series and at L=6 kg/m²s, G=0.3 kg/m²s, P=30 bar.

Figure 9.

Mean bulk temperature profiles for axial coordinate at different temperatures (L = 6 kg/m²s, G = 0.3 kg/m²s, P = 30 bar, d_p=2 mm)

As suggested from the simulations displayed in Fig. 9, the temperature increase predicted by the VOF multiphase model was always higher than the experimental value. As the reactor operation is now simulated at higher liquid flow rate, one should expect a reduction on the maximum temperature increase along the trickle-bed reactor. Indeed, the maximum temperature increase obtained at L=2 kg/m²s was 8.9% whereas it became 5.2% at L=6 kg/m²s meaning that a lower residence time did not allow for such heating degree as one observed for higher liquid phase residence time. At the lowest temperature, the averaged relative errors were -0.2, -0.3, -0.5, -0.5 and -0.6% for z=0.2, 0.4, 0.6, 0.8 and 1m, respectively. The magnitude of these decreased from the hotspot zone onwards with higher temperatures so we obtained -0.2, -0.1 and -0.2% for z=0.6, 0.8 and 1m.

4.4. Evaluation of reaction map fields

Prominent research on continuous catalytic wet oxidation of hazardous compounds has identified that hydrodynamics such as phase holdup, pressure drop, and liquid distribution are fundamental criteria for the efficient design of trickle beds. In this ambit and as long as the multiphase flow environment in trickle-bed reactors is strongly influenced by the local phase distribution and the velocity of the gas and the liquid contained in them, here we attempt to examine computationally the interstitial flow in environmentally-based trickle-bed reactors given the inherent difficulties of conventional measurement methods which mainly confer a global view on the trickle bed performance.

Two representative computational mappings of interstitial flow environment are depicted by Figs. 10 and 11 at L=6 kg/m²s, G=0.3 kg/m²s, and P=30 bar.

Figure 10.

Instantaneous snapshot of liquid phase distribution coloured by total organic carbon concentration (L=6 kg/m²s, G=0.3 kg/m²s, P=30 bar, T=200ºC, d_p=2 mm)

The former map corresponds to an instantaneous snapshot of liquid phase distribution coloured by total organic carbon concentration while the latter one shows an instantaneous snapshot of liquid phase distribution coloured by temperature at the same operating conditions.

The three-dimensional mappings allow us to ascertain the liquid flow consisting of rivulets and films, the former exhibiting poor liquid-solid contacting regions and the latter being good liquid-solid contacting. Indeed, this computational behaviour was experimentally observed by several authors who have systematized the rivulet/film classification in liquid distribution studies using computer-assisted tomography (Lutran et al., 1991; Ravindra et al., 1997).

Figure 11.

Instantaneous snapshot of liquid phase distribution coloured by liquid temperature (L=6 kg/m²s, G=0.3 kg/m²s, P=30 bar, T=200ºC, d_p=2 mm)

As one can observe, the liquid distribution is not uniform at the catalyst scale which identified computationally the so-called flow maldistribution of trickle beds. These strong axial and radial heterogeneities were observed experimentally by means of a magnetic resonance imaging technique to directly measure the flow in a pore space of a packed bed (Suekane et al., 2003) that in turn affected the spatial distribution of reaction parameters such as the ones illustrated in Figs. 10 and 11.

5. Conclusions

Advanced CFD simulation codes can be used to intensify and demonstrate new potentialities and shortcomings in the field of modern environmental detoxification technologies. In this work, we have developed an integrated experimental and theoretical framework to give further insights into consistent, safer and cleaner remediation technology of wastewaters contaminated with hazardous compounds by catalytic wet oxidation.

First, our case study encompassed the development of a CFD multiphase model which accounts for the chemical kinetics and interphase coupling hydrodynamics in the constitutive equations. On the optimization of mesh density and time step, the CFD model correctly handled the effect of these numerical solution parameters on the computations of TOC conversion and temperature. Second, particular attention was paid to the axial profiles of temperature and TOC depletion rates. The CFD simulations have shown promising results on how fluid dynamics can be correlated with chemical reaction, namely on the prediction of total organic carbon conversions functionalized by nominal liquid flow rates and attained at different temperatures. Finally, the interstitial flow was analysed by three-dimensional computational mappings of gas-liquid distribution as well as TOC and thermal representations highlighting the inhomogeneous reactive flow environment within the trickle-bed reactor. This methodology encompassed a novel integrated tool which empowered a secure estimation of reaction properties on the simulation of hazardous and complex multiphase flows with adjustable parameters.

Acknowledgments

The authors gratefully acknowledged the financial support of Fundação para a Ciência e Tecnologia, Portugal.