CHAPTER FIVE

PARTICLE-TO-FLUID MASS & HEAT TRANSFER

Convective transport at low & high pressure. Free convection effects in Supercritical Fluids

« Não há normas. Todos os homens são excepção a uma regra que não existe »

Fernando Pessoa (1888 – 1935) Image: Henri Cartier-Bresson, « Praying »

Chapter Five

One of the most important parameters needed in the design of packed bed systems is the particle-to-fluid transport coefficients. In the past six decades, a substantial amount of work has been devoted to the study of these parameters. Particle-to-fluid mass transfer studies were first carried out by Gamson et al. (1943), and Hurt (1943). They obtained mass transfer coefficients from measurements of the rates of evaporation of water from wet porous particles. Hurt (1943) also reported mass transfer coefficients derived from the measurement of rates of naphthalene sublimation. Since their pioneering work, a large number of experimental studies have been carried out on mass transfer coefficients in packed bed systems. Theoretical work has also been in progress. Pfeffer (1964) and Pfeffer and Happel (1964) applied a free surface cell model to the creeping flow region. LeClair and Hamielec (1968, 1970) proposed a zero vorticity cell model, and El-Kaissy and Homsy (1973) applied the free surface cell model, zero vorticity cell model and distorted cell model to a multiparticle system at low Reynolds number. Nishimura and Ishii (1980) also applied the free surface cell model to the study of mass transfer at high Reynolds numbers. These models which are based on different assumptions, generally give different and inconsistent values of particle-to-fluid mass transfer coefficients. Therefore, theoretical prediction of mass transfer coefficients is far from satisfactory. For the design and analysis of packed bed catalytic reactors, it is necessary to know the temperatures of the fluid and the catalyst particles in which the chemical reactions are taking place. In general, fluid temperature is measured with little difficulty, but the measurement of the solid surface temperature is not easy. This is particularly true of packed bed reactors. The particle temperature or temperature drop at the particle surface then has to be estimated in terms of the heat transfer coefficient between the particle and the fluid. Because of the importance of the particle-to-fluid heat transfer coefficient in packed bed reactors, a considerable effort has been made to evaluate this parameter. Experimental determinations of heat transfer coefficients for a wide variety of systems have been made using various experimental techniques, under either steady-state or unsteady-state conditions. An extensive review on experimental/theoretical works on particle-to-fluid mass and heat transfer can be found in Wakao and Kaguei (1982). Although there is broad information available of published experimental data for particle-tofluid forced convection heat and mass transfer at low pressure (see e.g. Wakao and Kaguei, 1982), scarce data are available for high pressure situations or supercritical packed bed reactors. Many authors have studied catalytic reactions using a supercritical fluid as a solvent (Poliakoff et al., 1996; Ramirez et al., 2004), but main interest has been centered in kinetics of chemical reactions under supercritical conditions. Other authors have studied mass transfer in packed beds under supercritical conditions applied to supercritical extraction (Debenedetti and Reid, 1986; Stüber et al., 1996). The recent study and development of new compact packed bed supercritical nuclear reactors opens a branch for the study of heat transfer phenomena under high pressure conditions (see e.g. Oka et al., 1994). In the following sections, analytical solutions of the steady-state temperature and composition profiles in a packed bed will be shown. In addition, a CFD simulation strategy for the estimation of forced convection mass and heat transfer coefficients at low pressure is discussed. Estimated transport parameters are compared to broadly accepted correlations. Mixed convection mass and heat transfer at high pressure is modeled and analyzed, and numerical results obtained are compared to previously published experimental data (Stüber et al., 1996), and a correlation for predicting heat transfer parameters in a supercritical packed bed reactor is presented (Guardo et al., 2006). The obtained numerical results validates the idea that the modified correlation presented by Guardo et al. (2006) can be used to describe heat transfer

- 93 -

Chapter Five

phenomena in a supercritical packed bed reactor under mixed convection regime at high pressures.

5.1.

Geometry, mesh design and CFD modeling

To properly design a mesh capable of capturing the transport mechanisms present in the study in detail, a dimensionless analysis of Navier-Stokes equations under simulation conditions was developed. The dimensionless equations corresponding to continuity, mass, momentum, and energy balances are detailed in Chapter two. The orders of magnitude of the dimensionless groups were estimated by taking physical-chemical property values for air from experimental data and empirical correlations available in the literature (Reid et al., 1987; Yaws, 1999; Poling et al., 2000). Reynolds number was calculated using particle diameter as characteristic length. Reynolds analogy was used to estimate values of Prt from Ret (White, 1991). For first analysis purposes, Prt was assumed as a constant value within the bed, justified in the fact that majority of experimental results shows a range of variations between 0.75 < Prt < 2 for air and water (Kays, 1994). Boundary (operating) conditions for each analyzed situation are shown in Table 5.1. Details on the dimensionless numbers used can be found in Appendix C. Results of the order of magnitude of the dimensionless groups are shown in Table 5.2.

Boundary condition

Low pressure

High pressure

Mass transfer simulations Circulating fluid C7H8 concentration at inlet, mol/m3 C7H8 concentration at particle surface (equilibria), mol/m3 Pressure, Pa Mass flow at the inlet, kg/m2· s Velocity at the inlet, m/s

CO2 0 5.95

120 – 190

101325 -4 7.5 x 10 – 7.5 x 10-1

9 – 9.2 x 106 0.015 – 0.100 -

Heat transfer simulations Circulating fluid Temperature at the inlet, K Temperature at the particle surface, K Pressure, Pa Mass flow at the inlet, kg/m2· s Velocity at the inlet, m/s Table 5.1.

Air 298 423 101325 -4 3 x 10 – 7.5 x 10-1

CO2 330 340 1 x 107 0.013 – 0.132 -

Boundary (operating) conditions for analyzed cases

In the case of particle-to-fluid transport at low pressure, dimensionless analysis allows identifying the problem as forced convection in laminar or turbulent flow. For the momentum balance it becomes clear that viscous forces decrease their contribution as Re increases. Inertial gravity forces increase their contribution as Re decreases, and pressure drop together with turbulent forces become the most important terms in the momentum balance at high Re. In energy and species balances, the convective and the diffusive term become important for the balances. For both balances, steady state analysis can be used.

- 94 -

Chapter Five

When analyzing a mass/heat transfer case at high pressure, the problem is identified as mixed (free + forced) convection in laminar or transitional flow. Body forces and pressure drop become the most important terms in the momentum balance. Turbulence forces contribution to all balance equations is negligible. In order to couple species/energy and momentum balance, unsteady state analysis is required.

Low Pressure

High Pressure

Re

100

101

102

103

100

101

Sr

102

101

100

10-1

105

104

Ma

10-6

10-5

10-4

10-3

10-6

10-5

Eu

103 - 104

101 - 103

10-1 - 101

10-1

106

104

Fr

10-4

10-4 - 10-2

10-2 - 100

100

10-10

10-8

Pr

10-1

10-1

10-1

10-1

100

100

Sc

100

100

100

100

100

100

Ec

10-10

10-10 - 10-8

10-8 - 10-6

10-6

10-16

10-14

ReT

101

100

100

100

101

100

PrT

10-1

10-1

10-1

10-1

10-1

10-1

ScT

103

102

101

100

103

102

Table 5.2.

Dimensionless groups’ magnitude orders for analyzed cases

In order to define the computational domain, a 44 spheres stacking with a sphere-to-tube diameters ratio of 3.923 was chosen for the geometrical model. Modeled geometry was constructed following the bottom-up technique (generating surfaces and volumes from nodes and edges) in order to control mesh size around critical points (i.e. particle-to-particle and particle-to-wall contact points). In this study, to include real contact points, the spheres were modeled overlapping by 0.5 % of their diameters with the adjacent surfaces in the geometric model. For further details on the geometrical model please refer to Guardo et al., (2004) and Chapter Four of this thesis.

5.2.

Model setup and analysis

Navier-Stokes equations together with species/energy balance were solved using commercially available finite volume code software Fluent® 5.x/6.x. For the low pressure models the fluid was taken to be incompressible, Newtonian, and in a laminar or turbulent flow regime. CO2, air and toluene at standard conditions were chosen as the simulation fluids. Incompressible ideal gas law for density, power law for viscosity and appropriated mixing rules were applied to the model for making these variables temperature and composition dependent. For the high pressure simulations, the fluid was taken to be Newtonian, in laminar flow regime and with variable density. CO2 and toluene were chosen in this case as the simulation fluids, and its

- 95 -

Chapter Five

properties at high pressure were incorporated to the solver code through User Defined Functions (UDF) and User Defined Equations (UDE). Peng-Robinson equation of state was used to calculate fluid’s density and heat capacity, Lucas method was used to calculate the viscosity, and thermal excess method was used to calculate the thermal conductivity (Reid et al., 1987). Appropriate mixing rules were applied to the model for making all variables compositiondependant. Details on the implementation of these UDE’s and UDF’s can be found in Appendix D. Under-relaxation factors for pressure, momentum and energy were initially set to 0.05, 0.1 and 0.2, respectively (Gunjal et al., 2005), and increased progressively after convergence until values of 0.2, 0.4 and 0.8, respectively. In the case of turbulent flow simulations, under-relaxation factors for turbulent quantities were set in 0.4. A first order discretization scheme for pressure, momentum and energy equations was used until convergence was achieved, and the results obtained were used as initial solution for a new simulation applying a second order discretization scheme for momentum and energy equations. Simulations were run in a Dell Precision 380 and Hewlett Packard Proliant DL385 workstations, and simulation times ranged from 1 to 240 hours depending on the studied case. Numerical convergence of the model was checked based on a suitable diminution of the normalized numerical residuals of all computed variables. For a more complete convergence checking the average static temperature or average species concentration at the bed outlet were also chosen as monitors depending on the case.

5.3.

Results and discussion

5.3.1.

Convective transport at low pressure

The main purpose of this set of simulations was to test CFD capabilities in a particle-to-fluid transport packed bed model. Standard correlations for Nu and Sh, and experimental data were selected as reference values to compare against the numerical results generated (Wakao et al., 1979, 1982). For each simulation, inlet velocity was varied (0.2 < Re < 1800), and the transport coefficients (Nu, Sh) were obtained. For the modeling in the turbulent flow region (Re > 300), previous experience shows that the Spalart-Allmaras turbulence model is a good choice (Guardo et al., 2005). In packed beds with constant temperature or species concentration at particle surface, the transport resistance resides only on the fluid side. If axial dispersion is taken into consideration, then the balance equations for steady-state analysis are as shown in Table 5.3. For each simulation, variable profiles along the bed were recorded and analyzed. With the collected numerical results, the transport coefficients (kc, h) were obtained from Equations [5.3-1] and [5.3-2]. From the obtained values for kc and h, Sh and Nu were computed and compared with the broadly accepted correlations and experimental data. Results for mass transfer are shown in Figure 5.1. Results for heat transfer are shown in Figure 5.2. In both mass and heat transfer simulation sets, four meshes were used. It can be noticed that for both mass and heat transfer situations, in the laminar and transition flow zone (Re < 300) the results do not show dependence on the mesh density. The refinement of the mesh is expressed in terms of Vcell / Vp ratio.

- 96 -

Chapter Five

Mass balance

Energy balance

2 dC 6k c (1 − ε ) (C f − C ps ) = α ax d C2 + u dx d p ⋅ε dx [5.3-1]

d 2T f 6h ⋅ (1 − ε ) (T f − T p ) = α ax 2 u + dx d p ⋅ ε ⋅ C p ⋅ ρ dx [5.3-2] dT f

Boundary conditions

u (C − C 0 ) = α ax

dC dx

u (T f − T0 ) = α ax

[5.3-1a]

at x = 0 (inlet)

[5.3-2a]

dx

at x = 0 (inlet)

dT f

dC =0 dx

[5.3-1b]

dx

at x = L (outlet) Table 5.3.

dT f

=0

[5.3-2b]

at x = xL (outlet)

Mass and energy balance equations applied to the packed bed model

At lower Reynolds numbers (Re < 10), mass and heat transfer results obtained shows that the fitting against Wakao et al. correlations (1979, 1982) is not good. This is particularly noticeable in the case of the heat transfer simulations (see Figure 5.2). For a single velocity condition, different meshes give results in a wide range of Nu and no relation with mesh density can be established in any case. Experimentally, one-shot measurements have demonstrated that no definite Nu values can be obtained at low Re (Wakao, 1976). The fact that any Nu value within the obtained ranges yields approximately the same value of αax suggests that at this low Re, particle-to-fluid heat transfer makes little contribution to the overall heat transfer in the system, which analytically has been reflected in an increased confidence range for the selected correlation at low Re (Shent et al., 1981). Using a laminar model, a good accuracy for Sh and Nu values was obtained within 10 < Re < 100, reinforcing the idea that a similar methodology can be applied to model mixed convection heat transfer at high pressure in the same Re range.

Sh

100

10

Sc = 1.3 1 1

10

100

1000

10000

Re VSerie10 cell / V p = 5.77 x 10

-4

= 5.74 x 10-4 Vvcell/vp cell / V p = 3.04 x 10

-5

Richardson & Szekely (1961)

Figure 5.1.

-4 VSerie11 cell / V p = 2.09 x 10

EXPERIMENTAL DATA (Reviewed by Experimental data

Wakao et al.)

VSerie12 cell / V p = 3.21 x 10

-5

Ranz & Marshall (1952)

Wakao & Kaguei (1982)

Sherwood number vs. Reynolds number for the low pressure convective mass transfer packed bed model

- 97 -

Chapter Five

100

Nu

RANGE OF EXPERIMENTAL DATA (Reviewed by Wakao et al.)

10

Pr = 0.9

1 0.1

1

10

Re -4

58933 V cell /Velements p = 5.77 x 10

-5 430051 = 3.04 x 10 V cell /V pelements

Figure 5.2.

104325 elements -4 V cell /V p = 2.09 x 10 Ranz & Marshall (1952)

100

1000

348065 elements -5 V cell /V p = 3.21 x 10 Wakao et al. (1979)

Nusselt number vs. Reynolds number for the low pressure convective mass transfer packed bed model

For higher values of Re there is a divergence between the results obtained for tested meshes in both heat and mass transfer simulations, probably due to the fact that at higher Re, turbulent transport term in the transport equation becomes more important. In RANS modeling, the balance equations possess a smooth exact solution, and the numerical solution approaches that solution as we refine the grid. The aim of grid refinement is numerical accuracy (Spalart, 2000; Shur et al., 2005). Therefore, for our specific study case, an accurate turbulence modeling requires a denser mesh around the particles surface in order to capture the involved turbulence phenomena and its related mixing enhancement in the boundary layer (Guardo et al., 2005). The results obtained with the finer meshes fit better the prediction of Wakao et al. (1979) for heat transfer, and Wakao and Kaguei (1982) for mass transfer in the turbulent flow zone (Re > 300), probably due to a better capture of the vorticity energetic scales associated effects.

5.3.2.

Convective transport at high pressure

A numerical modeling of a packed bed where mixed convection mass and heat transfer phenomena appear was also made. The same geometrical model as the aforementioned case was used but now the circulating fluid and the values of pressure and velocity were modified. The effects of density, flow rate and flow direction on mass and heat transfer are presented. Obtained numerical results for mass transfer were compared against the experimental results and the correlation presented by Stüber et al. (1996). Obtained numerical results for heat transfer were fitted, and a novel correlation based on the convective mass transfer correlation presented by Stüber et al. (1996) is proposed (Guardo et al., 2006a). An analogy between the obtained numerical data and previously reported heat transfer data (Guardo et al., 2006b) is also presented.

- 98 -

Chapter Five

Two simulation sets (upflow and downflow operation) were done, and for each simulation, inlet velocity was varied and Sh and Nu were obtained for 9 < Re < 100. A laminar flow model was used as viscous model for all simulations. Gravitational acceleration was added to the operating conditions of the model. Under-relaxation factors were set as previously explained for the aforementioned low pressure models. A first order discretization scheme for momentum and energy equations and a second order discretization scheme for pressure were used until convergence was achieved, and the results obtained were used as initial solution for a new simulation applying a second order discretization scheme for momentum and energy equations. Equations [5.3-1] and [5.3-2] were used to obtain the values of Sh and Nu for each studied case (see Table 5.3). Boundary conditions for the model are described in Table 5.1. The presence of high density gradients lead to the formation of hydrodynamical instabilities (i.e., stream differentiation, countercurrent flow, recirculation) caused by buoyancy effects (Benneker et al., 1998). In order to compute correctly the flow fields and the heat transfer phenomena in these zones, the optimal grid used for the prior model was locally refined for each case studied around the zones where high density gradients were found. Dimensionless analysis indicated that for the studied velocity range, laminar or transitional flow could be expected. For different mass and heat transfer simulations, local values of Gr, Pr, Sc and Re were calculated various points inside the boundary layer around the spheres in order to corroborate the transport mechanism present in the packed bed. Obtained values were compared with Metais-Eckert maps (Metais and Eckert, 1964) and are shown in Figure 5.3. As it can be observed, the main transport mechanism inside the packed bed under the different simulation conditions is the free convection in laminar/transitional flow regime.

1E+5

Forced Convection / Turbulent flow

Forced Convection / Laminar flow

1E+4

Free Convection / Turbulent flow

Re

1E+3

1E+2 Free Convection / Laminar flow

1E+1

1E+0 1E+2

1E+3

1E+4

1E+5

Gr · (Pr Figure 5.3.

or

1E+6

Transition: Laminar - Turbulent

1E+7

1E+8

1E+9

Sc) · (d/L)

Comparison between simulation data and Metais-Eckert maps

Churchill (1983) suggests that in a mixed convection heat transfer situation the transition between laminar and turbulent flow can be found at 1 × 109 < GrH < 1 × 1010. In order to corroborate this idea, the numerical results obtained for Nu for the different heat transfer situations analyzed were plotted in a flow map according to Churchill’s criteria (see Figure 5.4).

- 99 -

Chapter Five

It can be noticed that the results obtained follow the expected behavior according to the dimensionless analysis, and transition between laminar and turbulent flow shows up in the results by means of a high increase in Nu for alike values of GrH.

20 Laminar Flow

Transition

Turbulent Flow

Nu

15

10

5

0 1E+7

1E+8

1E+9

1E+10

1E+11

Gr H Assisting Flow Figure 5.4.

Opposing Flow

Flow map according to Churchill’s (1983) criteria

5.3.2.1. Effect of density gradients and flow stability Density gradients control the fluid movement at low velocities. Since the early experiments of Hill (1952), there has been a growth in the literature published on hydrodynamical instabilities due to density and viscosity variations in porous media (Homsy, 1987; Manickam and Homsy, 1995). Benneker et al. (1998) proved that the hydrodynamical instabilities caused by buoyancy effects in a packed bed when density increases with height lead to a significant increase in the flow axial dispersion within the bed. While simulating the high pressure data sets, the presence of hydrodynamical instabilities captured by the model caused numerical instability in the results. The presence of large density gradients may result in convergence problems when solving the momentum, continuity and pressure equations (Lakshminarayana, 1991). As mentioned before, in order to minimize the effects of flow instability on calculations, the mesh was locally refined in each case trying to avoid the presence of high density and velocity gradients within a single cell. Numerical instability increased as inlet mass flow imposed as boundary condition decreased. Figure 5.5 compares velocity fields (colored by density) for high and low pressure heat transfer situations. It can be seen that for the high pressure situation (Figure 5.5a) the presence of high density gradients (differences up to 40 kg/m3) lead to the differentiation of two streams: heavy (cold) fluid sinking into the bed, and light (hot) fluid floating out of the bed. This fact will be important when analyzing the effects of flow direction in free convection heat transfer. When analyzing the velocity–density profile for the low pressure situation (Figure 5.5b), it becomes

- 100 -

Chapter Five

clear that despite the fact a small density gradient is present in the results, flow field is controlled by the inlet velocity.

Figure 5.5.

Velocity field (scale colored by density, kg/m3) in an axial cut of the packed bed. (A) CO2 at 1 × 107 Pa and Re = 41; (B) Air at 101 325 Pa and Re = 531

- 101 -

Chapter Five

5.3.2.2. Effect of flow rate and flow direction As mentioned before, flow instability zones within the bed appear if a density gradient is present along the bed. The direction of this density gradient and the direction of the flow have strong effects on the axial dispersion and the mixing of the flow (Benneker et al., 1998). Either if the flow direction assists the transfer phenomena (when density decreases with height) or is opposed to it (when density increases with height), it can be noticed that, under certain conditions, flow direction have a direct effect on mass and heat transfer phenomena. Several authors have reported in experiments that under supercritical conditions, flow direction affects extraction rates increasing them when flow direction assists the mass transfer (e.g. Sovová et al., 1994; Stüber et al., 1996; Saiz et al., 2000; Germain et al., 2005). When gravity acceleration is added to the heat transfer model, the aforementioned phenomenology is captured by CFD analysis. Numerical runs were performed at various flow rates, constant pressure and temperatures. The fluid motion direction was also considered defining the fluid flowing in a direction opposite to gravity (upflow mode) or in its direction (downflow mode). When analyzing mass transfer, for the case of extraction of toluene from the particles surface, Figure 5.6 presents a solute mass fraction contour plot within the packed bed, and Figures 5.7 and 5.8 illustrates the initial mass extraction rates for different flow rates and flow directions. From the figures, two conclusions can be drawn. (1) Regardless of flow direction, the extraction rate increases strongly with flow rate (see Figure 5.7); and (2) Despite the strong dependence on flow rate, gravity-assisting flow enhances the extraction rate significantly (see Figures 5.6 and 5.8). This effect appears to be more pronounced at low Re. This behavior is certainly related to the presence of free convection and has been experimentally noted by other authors (Sovová et al., 1994; Stüber et al., 1996).

Figure 5.6.

Toluene mass fraction contour plot in an axial cut of the packed bed for [A] opposing flow and for [B] assisting flow with Re ≈ 20

- 102 -

Chapter Five

3.50E-05

C7H8 extracted [kg]

3.00E-05 2.50E-05 2.00E-05 1.50E-05 1.00E-05 5.00E-06 0.00E+00 0

30

60

90

120

150

Time [s] Re = 13 Figure 5.7.

Re = 20

Re = 36

Re = 78

Cumulative toluene extracted at the bed outlet vs. time for different Re in opposing flow regime

C7H8 extracted [kg]

2.50E-05 2.00E-05 1.50E-05 1.00E-05 5.00E-06 0.00E+00 0

30

60

90

120

150

Time [s] Opposing flow Figure 5.8.

Assisting flow

Cumulative toluene extracted at the bed outlet vs. time for different flow directions

- 103 -

Chapter Five

When analyzing heat transfer at high pressure, results obtained show a similar behavior than the aforementioned, reflecting that heat transfer coefficient values are significantly greater in assisting flow cases than those obtained for opposing flow cases (see figure 5.9). The heat transfer rising in assisting flow is favored by the temperature gradient direction (the same as flow direction), which increases flow mixing and heat transfer rates. The presence of an adverse density gradient (when temperature gradient direction is opposed to flow direction) diminishes heat transfer rates and lower coefficients are obtained.

20 16

Nu

12 8 4 0 40

50

60

70

80

90

100

Re Assisting flow Figure 5.9.

Opposing flow

Nusselt number vs. Reynolds number for assisting and opposing flow simulations

In a similar way, flow velocity has a directly proportional effect over heat transfer when forced convection takes place. An increase on flow velocity leads to an increase on kinetic energy. This fact generates a better mixing and a greater heat transfer coefficient (see Figure 5.2). When analyzing a mixed convection case, the obtained results indicate that despite the fact that the forced convection component of the phenomena follows this behavior, free convection component is almost constant for all cases (see Figure 5.10). Is clear that free convection component is intrinsically related to density and temperature fields and it is not affected by velocity. Therefore, the dominant dimensionless group describing this phenomenon will be GrH rather than Re. Independently of flow direction, for the covered range of Re the heat transfer coefficient value increases almost linearly with velocity (see Figure 5.9). This increase in the overall heat transfer is due exclusively to the increase of the contribution of forced convection to the total heat transfer rate.

- 104 -

Chapter Five

Nu

100

10

1 1

10

100

Re Free Convection

Figure 5.10.

Forced Convection

Nusselt number vs. Reynolds number for free and forced convection in assisting flow

5.3.2.3. Validating the numerically obtained mass transfer data Values of kc were evaluated from the mass balance equation (see Table 3), using numerical data and substance properties at corresponding operating conditions. From the obtained results for values of kc, Sh was computed and compared against experimental data and against the correlation presented by Stüber et al. (1996). The results of Sh vs. Re and the aforementioned comparison are given in Figure 5.11. As it can be seen, a good agreement between numerical results, experimental results and the correlation prediction is obtained.

Sh

100

10

1 10

CFD Figure 5.11.

100

Re Experimental

Stüber et al.

CFD obtained Sherwood number vs. Reynolds number for both upflow and downflow operations, and comparison against experimental data and correlations presented by Stüber et al. (1996).

- 105 -

Chapter Five

The contributions of free and forced convection to mass transfer could be evaluated following the indications of Stüber et al. (1996), finding also concordance between the numerically obtained results and the experimental results reported in that paper (see Figures 5.12 and 5.13).

Shfree/Sh (%)

100

10 10

Opposing flow

[(Sh - Sh0) ± (Shfree - Sh0)]/Sc

0.3

Figure 5.12.

Re Assisting flow

100

Experimental (Stüber et al.)

CFD obtained contribution of free convection to total mass transfer vs. Reynolds number and comparison against experimental data presented by Stüber et al. (1996).

100

10

1 10

100

Re Opposing flow Figure 5.13.

Assisting flow

Experimental (Stüber et al.)

CFD obtained contribution of forced convection to total mass transfer vs. Reynolds number and comparison against experimental data presented by Stüber et al. (1996).

- 106 -

Chapter Five

5.3.2.4. Correlating the numerically obtained heat transfer data Heat transfer coefficients at high pressure were obtained using a CFD simulation procedure similar to that applied for the mass transfer simulations presented in the prior sections of this chapter. The obtained numerical results validates the idea that a modified correlation based on the one proposed by Stüber et al. (1996) for mixed convection mass transfer under supercritical conditions in packed beds can be used to describe heat transfer phenomena in a supercritical packed bed reactor under mixed convection regime at high pressures. Such correlation can be useful when designing supercritical packed bed reactors, particularly in the design of directcycle supercritical-water-cooled fast nuclear reactors (Oka et al., 1994), nowadays under study and development process. Unfortunately, there is no experimental data available in order to validate the prior statement.

Nufree / Nu (%)

100

10 1

10

Re Assisting Flow

Figure 5.14.

100

Opposing Flow

CFD obtained contribution of free convection to total heat transfer as a function of Reynolds number

The aforementioned idea of a high pressure mixed convection heat transfer correlation can be noticed explicitly in Figure 5.14 and Figure 5.15. Figure 5.14 shows the contribution of free convection to the total Nusselt number as a function of the Reynolds number. As expected, low Re numbers correspond to transport dominated by natural convection, but for Re around 100, still 30% of transport is due to natural convection (for both assisting and opposing flow). Figure 5.15 shows a reduction of the numerical data points to a natural convection-free basis only depending on the Re numbers. While the contribution of free convection to total heat transfer rate decreases with Re (as its value remains almost a constant), forced convection increases, leading to greater heat transfer coefficients.

- 107 -

Chapter Five

[(Nu - Nu0 ) ± (Nufree - Nu0 )]/Pr

0.3

10.00

1.00

0.10 1

10

Assisting Flow Figure 5.15.

100

Re Opposing Flow

CFD obtained contribution of forced convection to total heat transfer as a function of Reynolds number

After validating the numerical results obtained against experimental data and empirical correlations (see Figures 5.11, 5.12 and 5.13), an analogy between mass transfer data presented in this work (CFD and experimental) and numerical heat transfer data reported (Guardo et al., 2006) can be done, establishing a relationship between the selected data. It can be noticed that with a similar simulation methodology heat and mass transfer coefficients at high pressure can be obtained, because the CFD code applies the same mathematical procedure to estimate transport of any scalar quantity (such as mass or energy). Results obtained for heat and mass transfer simulations show a good fitting when compared to the prediction of the aforementioned correlations for heat and mass transfer. The parity plot for predicted vs. obtained Nusselt and Sherwood numbers is given in Figure 5.16. In summary, the correlations recommended are:

Nu − Nu 0 = Nu forced m (Nu free − Nu 0 )

[5.3-3]

Taking Nu0 = 2 (Wakao and Kaguei, 1982), the correlations for free and forced convection are:

Nu free = Nu 0 + 0.001(GrH ⋅ Pr )

0.33

Pr 0.244

[5.3-3a]

Nu forced = 0.269 Re 0.88 Pr 0.3

[5.4-3b]

Valid for the following ranges of dimensionless numbers: 9 < Re < 96; 2.2 < Pr < 3.3; 1 × 108 < GrH · Pr < 4 ×1010.

- 108 -

Chapter Five

An advantage of this correlation is that a single set of parameters predicts heat transfer for assisting flow free convection (with the minus sign in Equation [5.3-3]) as well as opposing flow free convection (with plus sign in Equation [5.3-3]). Predicted values for Nusselt number show that estimation error is within that obtained with original mass transfer correlation, that showed an AARD = 18.9% (see Figure 5.17).

Coefficient from CFD

20

15

10

5

0 0

5

10

15

20

Coefficient from correlations Heat transfer Figure 5.16.

Mass transfer

Comparison of numerically obtained vs. predicted heat and mass transfer coefficients using the proposed correlations

25

20

NuT

15

10

5

0 1

2

3

4

5

6

7

8

9

10

11

Simulation run CFD

Figure 5.17.

Correlation [Eq. 5.3-3]

Estimated error for the high pressure mixed convection correlation against simulation runs

- 109 -

Chapter Five

CONCLUSIONS CFD proves to be a reliable tool when modeling convective mass and heat transfer phenomena in packed beds. It allows to analyze either free or forced convection situations and the obtained results can be compared qualitatively and quantitatively against previous published data. Forced convection at low pressure in a packed bed was simulated, and the influence of velocity over mass and heat transfer could be analyzed. For values of Re > 10 it was obtained a good agreement with the correlations presented by Wakao and Kaguei (1982) for mass transfer, and Wakao et al. (1979) for heat transfer. At lower Reynolds numbers (Re < 10), results shows that the fitting against correlation is not good. This situation is especially notorious when analyzing heat transfer. For a single velocity condition different meshes give results in a wide range of Nu and no relation with mesh density can be established in any case. No mesh sensitivity was noticed for the laminar and transition flow zones, but in the turbulent flow zone a good definition of the mesh around the particles surface is of primal importance in order to capture the turbulence vorticity energetic scales associated effects. Mixed (free+forced) convection at high pressure in a packed bed was also analyzed. For a supercritical fluid in laminar flow regime it was possible to study the effects of the density gradient, flow direction and velocity over mass and heat transfer. It was noticed that the presence of large density gradients conditioned the mesh influence over the numerical results when computing the mixing and the mass/heat transfer within the computed domain. An adverse density gradient generates hydrodynamical instabilities that produce an increase of the axial dispersion and a diminished mass/heat transfer rates. This fact caused numerical instability in the simulation process. In order to eliminate numerical instability, the mesh was locally refined in each case trying to avoid the presence of high density and velocity gradients along the bed. Influence of flow direction over mass/heat transfer was also analyzed and it was noticed that in assisting flow regime, greater heat transfer rates are obtained. The rising in the density along the bed height due to the adverse density gradient helps the axial dispersion to grow, obtaining less mixing within the bed and smaller heat transfer coefficients. Flow velocity also affects the heat transfer rate. The value of Sh and Nu increases almost linearly with flow velocity and this increase is due exclusively to the increase in the contribution of forced convection to the overall mass/heat transfer. Free convection is independent of flow velocity and its value remains almost constant within the studied velocity range. A novel correlation (Eq. [5.3-3]) is presented for estimating the free and forced convection effects and Nu from Re, Pr and GrH . Presented correlations are valid for 9 < Re < 96, 2.2 < Pr < 3.3, and 1 × 108 < GrH · Pr < 4 × 1010. The obtained numerical results validates the idea that the modified correlation presented (Guardo et al., 2006), based on the one proposed by Stüber et al. (1996) for mixed convection mass transfer under supercritical conditions in packed beds, can be used to describe heat transfer phenomena in a supercritical packed bed reactor under mixed convection regime at high pressures.

- 110 -

Chapter Five

REFERENCES Benneker, A., Kronberg, A., Westerterp, K., (1998). Influence of buoyancy forces on the flow of gases through packed beds at elevated pressures. AIChE Journal, 44, 263 – 270. Churchill, S., (1983). Single phase convective heat transfer. In: Schlünder, E. et al. (Eds.): Heat Exchanger Design Handbook. Hemisphere publishing corporation, New York, NY, USA. (Sections 2.5.7 to 2.5.10) Debenedetti, P., and Reid, R.C., (1986). Diffusion and mass transfer in supercritical fluids. AIChE Journal, 32, 2034 – 2046. El-Kaissy, M.M., and Homsy, G.M., (1973). A theoretical study of pressure drop and transport in packed beds at intermediate Reynolds numbers. Industrial & Engineering Chemistry Fundamentals, 12, 82 – 90. Gamson, B.W., Thodos, G., Hougen, O.A., (1943). Heat, mass and momentum transfer in the flow of gases through granular solids. Transactions of the American Institute of Chemical Engineers, 39, 1 – 35. Germain, J., del Valle, J., de la Fuente, J., (2005). Natural convection retards supercritical CO2 extraction of essential oils and lipids from vegetable substrates. Industrial & Engineering Chemistry Research, 44, 2879 – 2886. Guardo, A., Coussirat, M., Larrayoz, M.A., Recasens, F., Egusquiza, E., (2004). CFD Flow and heat transfer in nonregular packings for fixed bed equipment design. Industrial & Engineering Chemistry Research, 43, 7049 – 7056. Guardo, A., Coussirat, M., Larrayoz, M.A., Recasens, F., Egusquiza, E., (2005). Influence of the turbulence model in CFD modeling of wall-to-fluid heat transfer in packed beds. Chemical Engineering Science, 60, 1733 – 1742. Guardo, A., Coussirat, M., Recasens, F., Larrayoz, M.A., Escaler, X., (2006a). CFD study on particle-to-fluid heat transfer in fixed bed reactors: Convective heat transfer at low and high pressure. Chemical Engineering Science, 61, 4341 – 4353. Guardo, A., Coussirat, M., Recasens, F., Larrayoz, M.A., Escaler, X., (2006b). CFD studies on particle-to-fluid mass and heat transfer in packed beds: free convection effects in supercritical fluids. Chemical Engineering Science (submitted – in revision). Gunjal, P., Ranade, V., Chaudhari, R., (2005). Computational study of a single-phase flow in packed beds of spheres. AIChE Journal, 51, 365 – 378. Hill, S., (1952). Channeling in packed columns. Chemical Engineering Science, 1, 247 – 253. Homsy, G.M., (1987). Viscous fingering in porous media. Annual Reviews of Fluid Mechanics, 19, 271 – 311. Hurt, D.M., (1943). Principles of reactor design. Industrial & Engineering Chemistry, 35, 522 – 528. Kays, W.M., (1994). Turbulent Prandtl number. Where are we? Journal of Heat Transfer, 116, 284 – 295.

- 111 -

Chapter Five

Lakshminarayana, B., (1991). An assessment of computational fluid dynamic techniques in the analysis and design of turbomachinery. The 1990 Freeman scholar lecture. Journal of Fluids Engineering, 113, 315–352. LeClair, B.P., and Hamielec, A.E., (1968). Viscous flow through particle assemblages at intermediate Reynolds numbers. Steady-state solutions for flow through assemblages of spheres. Industrial & Engineering Chemistry Fundamentals, 7, 542 – 549. LeClair, B.P., and Hamielec, A.E., (1970). Viscous flow through particle assemblages at intermediate Reynolds numbers. Steady-state solutions for flow through assemblages of cylinders. Industrial & Engineering Chemistry Fundamentals, 9, 608 – 613. Manickam, O., and Homsy, G.M., (1995). Fingering instabilities in vertical miscible displacement flows through porous media. Journal of Fluid Mechanics, 288, 75 – 102. Nishimura, Y., and Ishii, T., (1980). An analysis of transport phenomena for multi-solid particle systems at higher Reynolds numbers by a standard Karman—Pohlhausen method—I Mass transfer. Chemical Engineering Science, 35, 1205 – 1209. Oka, Y., Koshizuka, S., Jevremovic, T., Okano, Y., (1994). Systems design of direct-cycle supercritical-water-cooled fast reactors. Nuclear Technology, 109, 1 – 10. Pfeffer, R., (1964). Heat and mass transport in multiparticle systems. Industrial & Engineering Chemistry Fundamentals, 3, 380 – 383. Pfeffer, R., and Happel, J., (1964). An analytical study of heat and mass transfer in multiparticle systems at low Reynolds numbers. AIChE Journal, 10, 605 – 611. Poliakoff, M., George, M.W., Howdle, S.M., (1996). Inorganic and related chemical reactions in supercritical fluids. In: Van Eldik, R. (Ed.), Chemistry under Extreme and Non-classical Conditions. Spektrum, Heidelberg, pp. 189–218 (Chapter 5). Poling, B., Prausnitz. J., O’Connell, J., (2000). The properties of gases and liquids. McGraw-Hill, New York. pp. 3.1 – 10.56. Ramírez, E., Recasens, F., Fernández, M., Larrayoz, M.A., (2004). Sunflower oil hydrogenation on Pd/C in SC propane in a continuous recycle reactor. AIChE Journal, 50, 1545 – 1555. Reid, R., Prausnitz, J., Poling, B., (1987). The properties of gases and liquids. McGraw-Hill, Boston. pp. 95 – 205. Saiz, S., Larrayoz, M.A., Trabelsi, F., Recasens, F., (2000). Procedimiento para la extracción de lanolina de la lana y planta para la realización de dicho procedimiento. Spain patent No. P200002461. Shent, J., Kaguei, S., Wakao, N., (1981). Measurements of particle-to-gas heat transfer coefficients from one-shot thermal responses in packed beds. Chemical Engineering Science, 36, 1283 – 1286. Shur, M., Spalart, P., Squires, K., Strelets, M., Travin, A., (2005). Three dimensionality in Reynolds-averaged Navier–Stokes solutions around two-dimensional geometries. AIAA Journal, 43, 1230 – 1242. Sovová, H., Kučera, J., Jež, J., (1994). Rate of the vegetable oil extraction with supercritical CO2 – II. Extraction of grape oil. Chemical Engineering Science, 49, 415 – 420.

- 112 -

Chapter Five

Spalart, P.R., (2000). Strategies for turbulence modelling and simulations. International Journal of Heat and Fluid Flow, 21, 252 – 263. Stüber, F., Vázquez, A.M., Larrayoz, M.A., Recasens, F., (1996). Supercritical fluid extraction of packed beds: external mass transfer in upflow and downflow operation. Industrial & Engineering Chemistry Research, 35, 3618–3628. Wakao, N., (1976). Particle-to-fluid transfer coefficients and fluid diffusivities at low flow rate in packed beds. Chemical Engineering Science, 31, 1115 – 1122. Wakao, N., Kaguei, S., Funazkri, T., (1979). Effect of fluid dispersion coefficients on particle-tofluid heat transfer coefficients in packed beds: correlation of Nusselt numbers. Chemical Engineering Science, 34, 325 – 336. Wakao, N., and Kaguei, S., (1982). Heat and mass transfer in packed beds. Gordon and Breach Science Publishers, New York. pp. 138 – 160; 264 – 295. White, F., (1991). Viscous fluid flow. McGraw-Hill, New York. pp. 482 – 542. Yaws, C., (1999). Chemical properties handbook. McGraw-Hill, New York. pp. 1 – 55.

- 113 -

CHAPTER SIX

HIGH–PRESSURE HETEROGENEOUS REACTIONS

CFD modeling on the supercritical hydrogenation of sunflower oil

« Verdrink niet in je eigen inspiratie en verbeelding; wordt geen slaaf van je eigen model »

Vincent van Gogh, « Het zonnebloemen » (1888)

Chapter Six

Partial hydrogenation of fats and oils is an important process in the food industry because not only increases oxidative stability of the final product when compared to raw materials, but also changes the physical characteristics of it, according to their final application (e.g. margarines or shortenings). In the conventional process, low reaction rates and the formation of undesirable by-products (such as trans fatty acids content about 40 wt%) are consequences of the low solubility of H2 in the oil and the high mass transfer resistance for the hydrogen in the liquid phase (Farrauto and Bartholomew, 1997). Supercritical technology has proven to be a reliable alternative to the conventional hydrogenation process (Härrod et al., 2001; Tacke et al., 2003) resulting increasing the rate of reaction and reduced trans isomer levels. It also provides a clean, economic and environmental friendly process. Porous solid catalysts used for gas catalytic reactions have specific surface areas of tens to hundreds of square meters per gram. This enormous amount of surface area results mainly from the fine interconnecting pores in the catalyst pellets. If a chemical reaction is very fast, it proceeds at the external surface of the pellet. If, however, the reaction is very slow, the reactant gas may diffuse deep into the pores of pellet, even to the center of the pellet, and the chemical reaction takes place everywhere uniformly in the pellet. In the laboratory, reaction rate determined directly by measurements using a differential reactor, for example, is the overall rate. The overall rate constant does not necessarily mean the intrinsic chemical reaction rate constant. For the design of industrial packed bed reactors, one needs to know the overall reaction rate, not the intrinsic chemical reaction rate. The overall rate is governed not only by the chemical reaction, but also by the diffusion rate through the pores inside the catalyst pellet as well as at the pellet’s external surface. If we simply measure activation energy from the overall reaction rate constants, the activation energy may differ from that of the intrinsic chemical reaction. The importance of diffusion is often underestimated by some catalyst chemists. Wheeler (1951, 1955) made an extensive study on the role of pore diffusion in catalysis. Also, Petersen (1965), Smith (1970), Satterfield (1970), Jackson (1977), and Dullien (1979) have reviewed the subject of pore diffusion associated with chemical reaction well. The review article of Youngquist (1970) will also help readers understand the basic principles of diffusion and reaction in a porous catalyst. In previous work done by our research group (Ramírez et al., 2004), an experimental design methodology and a response surface methodology were used to achieve optimum hydrogenation conditions for Pd catalyst in supercritical propane in a continuous well-mixed reactor. The results showed that it is possible to obtain a hydrogenated fat with 2 – 3 wt% trans C18:1 content and an iodine value (IV) of about 70 (starting from an initial value of 130). External resistance to mass transfer was made negligible using a high flow velocity over the catalyst bed (Re ≈ 6000). Furthermore, studies applying the Hashimoto et al. (1971) scheme provided information on the kinetics for the supercritical hydrogenation on a commercial-size Pd/C catalyst. In later studies (Ramírez et al., 2006), the intra-particle diffusion mechanisms of the triglycerides and hydrogen in the hydrogenation were studied using supercritical propane for 2 % Pd on activated carbon. The method of measuring the effective diffusion coefficients was first to obtain the intrinsic kinetic constants on small-size catalyst particles (≈ 0.5 mm) in the absence of diffusional limitation and then using the intrinsic kinetic constants, to derive the diffusivities in diffusion-limited reaction runs using larger catalyst sizes (up to 2 mm). The results showed that while hydrogen is transported by bulk pore diffusion, the oil seem to diffuse by surface diffusion. Diffusivity for H2 is about 10 times greater than that for triglycerides.

- 117 -

Chapter Six

The aim of the work presented in this chapter was to apply Computational Fluid Dynamics (CFD) to the study of the catalytic hydrogenation of sunflower oil in the presence of a supercritical solvent. A 2D CFD model of a single Pd-based catalyst pellet is presented. Intraparticle concentration profiles for all species present in the mixture (oil triglycerides and hydrogen) are obtained and compared against experimental results. Different particle sizes are studied and external mass transfer and intra-particle diffusional effects are analyzed. Verified kinetic constants and intra-particle diffusion coefficients were fed into a 3D packed bed reactor model, and conversion profiles are obtained.

6.1.

Geometry, mesh design and CFD modeling

To properly design a mesh capable of capturing the transport mechanisms present in the study in detail, a dimensionless analysis of Navier-Stokes equations under simulation conditions was developed. The dimensionless equations corresponding to continuity, mass, momentum, and energy balances are detailed in Chapter Two. The orders of magnitude of the dimensionless groups were estimated by taking physical-chemical property values for air from experimental data and empirical correlations available in the literature (Reid et al., 1987; Yaws, 1999; Poling et al., 2000). Reynolds number was calculated using particle diameter as characteristic length. Reynolds analogy was used to estimate values of Sct and Prt from Ret (White, 1991). Boundary (operating) conditions for each analyzed situation are shown in Table 6.1. Values for operating conditions were taken from experimental data published by Ramírez et al., (2006). Details on the dimensionless numbers used can be found in Appendix C. Results of the order of magnitude of the dimensionless groups are shown in Table 6.2.

Boundary condition

Single pellet model

Circulating fluid

Propane + H2 + Sunflower oil

Mixture proportions (wt %)

95 : 4 : 1

Temperature, K

457 – 473

Pressure, Pa

2 x 107

Flow model

Laminar; κ − ε

Laminar

0.13 – 1.32

1 x 10-3 – 2 x 10-2

Velocity at the inlet, m/s Particle diameter, mm L inlet concentration, mol/m3 O inlet concentration, mol/m3 E inlet concentration, mol/m3 S inlet concentration, mol/m3 H2 inlet concentration, mol/m3 Table 6.1.

Packed bed model

0.10 0.47

0.92

2.00

0.47

0.179669

0.217710 0.171589

37.73

0.075731

0.094676 0.057084

13.32

0.013177

0.018240 0.023614

0.00

0.108529

0.056489 0.120502

2.57

1.443499

1.5590

1.450727

Boundary (operating) conditions for analyzed cases

- 118 -

242.28

Chapter Six

Magnitude order

Table 6.2.

Re

100

101

102

103

104

Sr

10-1

10-2

10-3

10-4

10-3

Ma

10-6

10-5

10-4

10-3

10-3

Eu

1010

108

106

104

104

Fr

10-2

10-1

100

101

100

Pr

100

100

100

100

100

Sc

10-1

10-1

10-1

10-1

10-1

Ec

10-10

10-8

10-6

10-4

10-4

ReT

10-3

10-2

10-1

100

101

PrT

103

103

103

103

103

ScT

102

102

102

102

102

Dimensionless groups’ magnitude orders for analyzed cases

In the case of particle-to-fluid transport, dimensionless analysis allows identifying the problem as forced convection in laminar or turbulent flow. For the momentum balance it becomes clear that viscous forces decrease their contribution as Re increases. Inertial gravity forces increase their contribution as Re decreases, and pressure drop together with turbulent forces become the most important terms in the momentum balance at high Re. In energy and species balances, the convective and the diffusive term become important for the balances. For both balances, steady state analysis can be used. In order to define the computational domains, 2 geometrical models were created. A 2D axisymmetric simplification of a single spherical catalyst pellet (see Figure 6.1) was built to study mass transfer and intra-particle diffusional effects. 4 different particle sizes were selected for this study (Dp = 0.1, 0.47, 0.9205 and 2 mm). A simple cubic sphere stacking (see Figure 6.2) was selected as the base unit cell in order to simulate a packed bed reactor. For the simulations the fluid was taken to be Newtonian, isothermal (457 - 473 K) and in a laminar (10 < Rep < 100 for the packed bed model) or turbulent flow regime (u = 1.3 m/s for the single sphere model). A mixture of propane, hydrogen and sunflower oil (95 : 4 : 1) was chosen as the simulation fluid.

Figure 6.1. Mesh detail for the 2D single pellet geometrical model

- 119 -

Chapter Six

Figure 6.2.

Mesh detail for the 3D packed bed unit cell geometrical model

6.2.

Model analysis and setup

CFD has proven to be a reliable tool when modeling chemical reactors. Several reviews on the applicability of CFD methodology to chemical reactors design have been published in the last years (Bode, 1994; Harris et al., 1996; Kuipers and van Swaaij, 1998; Ranade, 2002). In this chapter, in order to study the transport and reaction mechanisms present in the supercritical hydrogenation of sunflower oil, Navier-Stokes equations together with the κ - ε turbulence model (when necessary) and a convection/diffusion mass balance were solved using a CFD commercially available finite element code software, Comsol Multiphysics 3.x. Simulations were run under isothermal (457 – 473 K) conditions. Steady state analysis was chosen for the simulations, following the guidelines obtained with the dimensionless analysis. The kinetic and diffusional models proposed by our research group (Ramírez et al., 2004, 2006) for the supercritical hydrogenation of sunflower oil were also incorporated within the equations set of the CFD solver. A schematic representation of the kinetic model used is shown in Figure 6.3. Kinetic constants used for simulations can be seen in Table 6.3. Values of molecular and effective intra-particle diffusivities applied to the CFD model can be seen in Table 6.4. For further insight in the kinetic and diffusional model, please refer to Ramírez (2005).

Figure 6.3.

Schematic representation of the kinetic model used for the supercritical hydrogenation of sunflower oil

- 120 -

Chapter Six

Temperature, K

Parameter

457

473

0.071

0.1

1 x 10-6

1 x 10-6

0.352

0.426

0.603

0.933

6.77 x 10-3

0.116

7.13 x 10-3

0.02

k12 mol-0.5· m4.5· (kg Pd)-1· s-1

k13 mol-0.5· m4.5· (kg Pd)-1· s-1

k23 mol-0.5·

m4.5·

mol-0.5·

m4.5·

(kg Pd)-1· s-1

k32 (kg Pd)-1· s-1

k24 mol-1· m6· (kg Pd)-1· s-1

k34 mol-1· m6· (kg Pd)-1· s-1 Table 6.3.

457 K Oil

H2

473 K Oil

H2

Table 6.4.

Fitted parameters values for the kinetic model

Particle diameter, mm Molecular diffusivity, m2/s Effective diffusivity, m2/s Molecular diffusivity, m2/s Effective diffusivity, m2/s

Particle diameter, mm Molecular diffusivity, m2/s Effective diffusivity, m2/s Molecular diffusivity, m2/s Effective diffusivity, m2/s

0.47

0.9205

2

2.94 x 10-8

2.94 x 10-8

2.94 x 10-8

2.94 x 10-8

1.76 x 10-8

8.82 x 10-9

6.4 x 10-8

6.4 x 10-8

6.4 x 10-8

6.4 x 10-8

3.84 x 10-8

1.92 x 10-8

0.47

0.9205

2

9.8 x 10-8

9.8 x 10-8

9.8 x 10-8

9.8 x 10-8

5.88 x 10-8

2.94 x 10-8

3 x 10-7

3 x 10-7

3 x 10-7

3 x 10-7

1.8 x 10-7

9 x 10-8

Fitted molecular and effective diffusion coefficients for hydrogenation species on 2% Pd/C catalyst (Dp range = 0.47 - 2 mm)

Simulations were run in a Dell Precision 380 workstation, and simulation times ranged from 1 to 4 hours depending on the studied case. Numerical convergence of the model was checked based on a suitable diminution of the normalized numerical residuals of all computed variables. Concentration contour fields (see Figure 6.4) for all species present in the model were obtained and the numerical data recorded in order to analyze conversion, reaction velocities and mass transfer.

- 121 -

Chapter Six

Figure 6.4.

(A) Concentration contour plot for Linoleic fatty acid (C18:2) in a single 2 % Pd catalyst pellet model (Dp = 0.9205 mm); (B) Concentration contour plots for fatty acids in a 3D packed bed model (Re = 100). Concentrations expressed in mol/m3

- 122 -

Chapter Six

6.3.

RESULTS AND DISCUSSION

6.3.1.

Validation of numerical results

In order to validate the numerical results obtained, a set of 2D simulations of particle-to-fluid mass transfer and reaction in a single catalyst pellet (emulating the experimental conditions proposed by Ramírez et al. (2006)) was made, and the results obtained were compared against the experimental results aforementioned. Intra-particle concentration profiles and intra-particle reaction rates were obtained for Dp = 0.47, 0.9205 and 2 mm, for which validation was possible due to the availability of experimental data.

0.24

1.6

0.18

1.2

0.12

0.8

0.06

0.4

CH2, surface [mol/m3]

Coil, surface [mol/m3]

First of all, it was necessary to verify the absence of gradients in the catalyst surface, since the experimental data was obtained in a Robinson-Mahoney gradientless reactor (Ramírez et al., 2004, 2006). For all particle diameters studied, concentrations (see Figure 6.5) and reaction velocity profiles (see Figure 6.6) along the surface arc were recorded. In all cases analyzed, flat variables profiles were obtained, assuring the absence of gradients along the catalyst surface.

Dp = 0.47 mm 0 -1.00

-0.50

0.00

0 1.00

0.50

Surface arc length [dimensionless] Linoleic Figure 6.5.

Oleic

Elaidic

Stearic

Hydrogen

Verification of the gradientless condition at the catalyst surface. Species concentration on catalyst surface for Dp = 0.47 mm

Once verified the aforementioned, intra-particle concentration profiles and intra-particle reaction rates were obtained and compared against experimental data previously published by our research group (Ramírez et al., 2006), obtaining for all cases a good agreement between the numerical and the experimental data (see Figures 6.7 and 6.8). Magnitude of relative errors was evaluated for the developed simulations, finding that this value was within the magnitude order of the experimental error previously reported by our research group. A trial for minimizing the relative error by means of optimizing the kinetic constants of the model and the effective diffusivities was also made, but an appreciable change was not observed in the concentration profiles obtained, remaining the magnitude order of the relative errors similar to those obtained when the original kinetic and diffusional model was used (Guardo et al., 2006).

- 123 -

0.03

-0.08

0.015

-0.11

0

-0.14

-0.015

-0.17

rH2, surface [mol/m3·s]

roil, surface [mol/m3·s]

Chapter Six

Dp = 2 mm -0.03

-0.2 -1

-0.5

0

0.5

1

Surface arc length [dimensionless] Linoleic Figure 6.6.

Oleic

Elaidic

Stearic

Hydrogen

Verification of the gradientless condition at the catalyst surface. Reaction rates on catalyst surface for Dp = 2 mm

1.6

0.25

1.4

0.8 0.1

0.6 0.4

0.05

3

0.2

Dp = 0.47 mm 0

0 -1

-0.5

0

Particle radius Linoleic Linoleic - EXP Hydrogen

Figure 6.7.

CH 2

Coil

1

0.15

[mol/m ]

1.2

3

[mol/m ]

0.2

0.5

1

[dimensionless]

Oleic Oleic - EXP Hydrogen - EXP

Elaidic Elaidic - EXP

Stearic Stearic - EXP

Validation of numerical data obtained. Intra-particle species concentration profile for Dp = 0.47 mm. Experimental data taken from Ramírez et al. (2006)

- 124 -

Chapter Six

0.0045

-0.005 -0.01

3

[mol/m ·s]

-0.015

0.0015 0

-0.025

-0.0015

r H2

-0.02

roil

3

[mol/m ·s]

0.003

-0.03

Dp = 2 mm -0.003

-0.035 -1

-0.5

0

Particle radius Linoleic Linoleic - EXP Hydrogen

1

[dimensionless]

Oleic Oleic - EXP Hydrogen - EXP

Elaidic Elaidic - EXP

Stearic Stearic - EXP

Validation of numerical data obtained. Intra-particle reaction rates for Dp = 2 mm. Experimental data taken from Ramírez et al. (2006)

0.25

1.5

0.2

1.4

0.15

1.3

0.1

1.2

0.05

1.1

CH 2

Coil

Figure 6.8.

0.5

1

0 0

0.2

0.4

0.6

0.8

1

Particle radius [dimensionless] Linoleic 0.1mm Linoleic 0.47mm Hydrogen 0.1mm

Figure 6.9.

Oleic 0.1mm Oleic 0.47mm Hydrogen 0.47mm

Elaidic 0.1mm Elaidic 0.47mm

Stearic 0.1mm Stearic 0.47mm

Comparison of the intra-particle concentration profiles for Dp = 0.1 mm (numerical) and Dp = 0.47 mm (experimental – Ramírez et al., 2006) in a single 2 % Pd catalyst pellet

- 125 -

Chapter Six

6.3.2.

Verification of catalyst effectiveness at low particle sizes

Due to the difficulty for obtaining experimental data at very low particle sizes, a CFD simulation extrapolating the catalyst evaluation for Dp = 0.1 mm was also made in order to corroborate the catalyst effectiveness. It can be seen in Figure 6.9 that for a particle with Dp = 0.1 mm, numerically obtained concentration profiles were almost identical to those reported experimentally for Dp = 0.47 mm (Ramírez et al., 2006), reassuring the idea that for a particle diameter of 0.47 mm or lower, the effectiveness factor of the catalyst is close to the unity.

6.3.3.

Particle-to-fluid mass transfer coefficients estimation

For each simulation, species concentration contour plots (Figure 6.4) were analyzed, and mass flux through the particle surface was recorded. With this data the surface local mass transfer coefficients (k) could be determined:

J i = k i ⋅ Ae ⋅ (C i , surface − C i ,∞ )

[6.3-1]

Figure 6.10 shows the surface local mass transfer coefficients for oil and Figure 6.11 shows the surface local mass transfer for H2 obtained by means of the CFD simulations. It can be seen that the obtained mass transfer coefficient values are almost constant in the upstream surface sector, while in the downstream surface sector the mass transfer coefficient value is strongly diminished. This fact must be related with the formation of an stagnation/recirculation flow point downstream the catalyst pellet, affecting specially the mass fluxes over the catalyst surface. 0.05

0.10 mm

[m/s]

0.03

koil, local

0.04

0.02

0.47 mm

0.9205 mm

0.01

2 mm 0

-1.00

-0.50

0.00

0.50

1.00

Surface arc length [dimensionless] Figure 6.10.

Surface local mass transfer coefficients for oil fatty acids for the supercritical hydrogenation of sunflower oil at 20 MPa and 473 K

- 126 -

Chapter Six

0.09

0.10 mm

kH2, local [m/s]

0.07

0.05

0.47 mm 0.04

0.9205 mm

0.02

2 mm 0.00

-1.00

-0.50

0.00

0.50

1.00

Surface arc length [dimensionless] Figure 6.11.

Surface local mass transfer coefficients for hydrogen for the supercritical hydrogenation of sunflower oil at 20 MPa and 473 K

From the values of k, the Sherwood number (Sh) was computed (see Table 6.5) and compared with the theoretical solution proposed by Ranz and Marshall (1952) for particle-to-fluid mass transfer in a single sphere (gas at low pressure), finding that obtained values for the Sherwood number are higher than those predicted by the aforementioned theoretical correlation (see Figure 6.12). When comparing the obtained results with an extrapolation for the correlation proposed by Tan et al. (1988) for supercritical extraction in packed beds or with an extrapolation of the experimental results obtained by Brunner (1985) for a high pressure bubble column, it can be noticed that the obtained values are lower than the extrapolation predictions. The absence of experimental data at such high Reynolds numbers difficult the validation tasks of the obtained numerical data. Nevertheless, the obtained numerical data suggests that the particle-to-fluid mass transfer coefficients for the supercritical hydrogenation of sunflower oil lie within the values for a gas (Ranz and Marshall, 1952) and a liquid (Brunner, 1985), which is coherent with the definition of a supercritical fluid (Guardo et al., 2006).

Dp, mm 0.10 0.47 0.9205 2.00 Table 6.5.

Re

L

Sherwood number O E S

H

647 44.16 44.16 44.18 44.18 27.61 3044 113.43 113.47 113.56 113.78 64.11 5962 208.14 208.32 207.11 209.28 119.85 12954 628.94 621.55 622.79 620.35 360.41

Computed values for Sherwood number for oil fatty acids and hydrogen in the supercritical hydrogenation of sunflower oil at 20 MPa and 473 K at different particle sizes

- 127 -

Chapter Six

100000 Brunner (1985)

10000

Sh

1000

THIS WORK

Tan et al. (1988) 100 10 Ranz & Marshall (1952) 1 10

100

1000

10000

Re p HYDROGEN

100000

OIL

Figure 6.12.

Sherwood number vs. Reynolds number for different particle diameters on the supercritical hydrogenation of sunflower oil at 20 MPa and 473 K

6.3.4.

Temperature effects on external mass transfer coefficients

In order to evaluate the temperature effects on the external mass transfer coefficients for the supercritical hydrogenation of sunflower oil, a simulation set at two different temperatures (457 and 473 K) was done. Particle size of Dp = 0.47 mm was chosen for the simulations because of being considered the optimal particle size for the hydrogenation reaction, as proven with the diffusional model and the intra-particle concentration profiles described in prior sections of this chapter. Values for the kinetic parameters and the values for the diffusion coefficients were set as show in Table 6.3 (for the kinetic model) and Table 6.4 (for the diffusional model). From the numerical results obtained, surface local mass transfer coefficients (see Figure 6.13) and Sherwood numbers (see Table 6.6) for the temperatures selected were estimated, following a procedure similar to that explained priory on this chapter. It can be noticed that a temperature increase is reflected in a direct increase of the surface local mass transfer coefficients. When expressing these numerically obtained mass transfer coefficients in a dimensionless way, it can also be noticed that mass transfer coefficients increase with Reynolds.

Temperature, K 457 473 Table 6.6.

Re

L

Sherwood number O E S

H

3972 202.41 202.24 198.86 196.80 139.30 3044 113.43 113.47 113.56 113.78 64.11

Computed values for Sherwood number for oil fatty acids and hydrogen in the supercritical hydrogenation of sunflower oil at different temperatures (particle size, Dp = 0.47 mm)

- 128 -

Chapter Six

0.03

0.045

473 K 0.03375

H2 OIL

0.015

0.0225

457 K

H2

0.0075

0.01125

OIL

kH2, local [m/s]

koil, local [m/s]

0.0225

Dp = 0.47 mm 0 -1.00

-0.50

0.00

0.50

0 1.00

Surface arc length [dimensionless] Figure 6.13.

Temperature effects over the surface local mass transfer coefficients for the supercritical hydrogenation of sunflower oil

6.3.5. Superficial velocity effects on external mass transfer coefficients In order to evaluate the superficial velocity effects on the external mass transfer coefficients for the supercritical hydrogenation of sunflower oil, a simulation set at six different inlet flow velocities (u = 2.65, 1.32, 0.66, 0.132, 0.066 and 0.0132 m/s) was done. Particle size of Dp = 0.47 mm was chosen for the simulations because of being considered the optimal particle size for the hydrogenation reaction, as explained in prior sections of this chapter. Values for the kinetic parameters and the values for the diffusion coefficients were set as show in Table 6.3 (for the kinetic model) and Table 6.4 (for the diffusional model).

Superficial velocity, m/s 0.0132715 0.0663575 0.1327150 0.6635750 1.3271500 2.6500000 Table 6.7.

Re

L

Sherwood number O E S

30 8.75 8.75 8.75 8.75 152 19.54 19.54 19.54 19.54 304 30.40 30.24 30.43 30.33 1522 111.76 111.61 111.78 111.69 3044 113.43 113.47 113.56 113.78 6077 114.63 114.43 114.68 114.54

H 5.73 12.72 17.54 63.89 64.11 64.89

Computed values for Sherwood number for oil fatty acids and hydrogen in the supercritical hydrogenation of sunflower oil at different inlet flow velocities (particle size, Dp = 0.47 mm)

- 129 -

Chapter Six

From the numerical results obtained, surface local mass transfer coefficients (see Figures 6.14 and 6.15) and Sherwood numbers (see Table 6.7) for the temperatures selected were estimated, following a procedure similar to that explained priory on this chapter.

kH2, local [m/s]

0.045

0.03

0.015

0 -1

-0.5

0

0.5

1

Surface arc length [dimensionless] v = 2.65 m/s v = 0.13 m/s Figure 6.14.

v = 1.32 m/s v = 0.066 m/s

v = 0.66 m/s v = 0.013 m/s

Surface velocity effects on the hydrogen local mass transfer coefficients for the supercritical hydrogenation of sunflower oil.

0.025

koil, local [m/s]

0.02 0.015 0.01 0.005 0 -1

-0.5

v = 2.65 m/s v = 0,13 m/s Figure 6.15.

0

0.5

Surface arc length [dimensionless] v = 1,32 m/s v = 0.066 m/s

1

v = 0,66 m/s v = 0.013 m/s

Surface velocity effects on the oil fatty acids local mass transfer coefficients for the supercritical hydrogenation of sunflower oil

- 130 -

Chapter Six

From Figures 6.14 and 6.15 and the numerical data shown in Table 6.7 it can be concluded that for high superficial velocities (u > 0.66 m/s), there is no significant effect of the superficial velocity over mass transfer phenomena. This behavior is in complete agree with the expected behavior for a well-stirred reactor. A similar situation has been studied experimentally by Ramírez et al., (2004) and Ramírez (2005), who reported for a CSTR supercritical hydrogenation reactor that an increase in operating stirrer speed (and in consequence, an increase in superficial velocity) from 52 rad/s to 262 rad/s produces no significant increase in the reaction rates of the hydrogenation. When taking a look at the surface local mass transfer coefficients (shown in Figures 6.14, 6.15 and Table 6.7) at low superficial velocities (similar velocities to those used in supercritical packed bed reactors) it can be noticed that there is a clear effect of velocity over mass transfer. Obtained mass transfer coefficients are significantly lower than those obtained for a well-stirred reactor. It was noticed that surface reaction velocities decreased with surface velocity, and that the total mass flux over particle surface increased with the surface velocity. The diminution of the surface reaction rates and the increase of the surface mass fluxes can significantly affect the surface mass transfer coefficients.

6.3.6. Correlating the numerically obtained mass transfer data Tan et al., (1988) correlated experimentally obtained mass transfer data for a fluid-solid system in a supercritical extractor, fitting their data with an equation of the type:

Sh = a ⋅ Re x ⋅ Sc y

[6.3-2]

Following the aforementioned model, numerically obtained Sherwood numbers (shown in Tables 6.5, 6.6 and 6.7) were correlated as a function of Reynolds and Schmidt numbers, in order to present a simple set of equations capable to predict mass transfer coefficients in the supercritical hydrogenation of sunflower oil. Figure 6.16 shows the obtained correlations for the numerically obtained mass transfer data. In summary, suggested correlations for predicting the oil fatty acids and hydrogen mass transfer coefficients in the supercritical hydrogenation of sunflower oil are as follows:

Shoil = 0.5192 ⋅ Re 0.69 ⋅ Sc 0.33

[6.3-3]

ShH 2 = 0.4796 ⋅ Re 0.6765 ⋅ Sc 0.33

[6.3-4]

Valid for the following range of dimensionless numbers: 30 < Re < 12954 0.7 < Sc < 2.1 AARD = 2.8 %

- 131 -

Chapter Six

1000

Sh oil = 0.5192 · Re

· Sc

0.33

[6.3-3]

Sh

100

0.69

10 Sh H 2 = 0.4796 · Re

0.6765

· Sc

0.33

[6.3-4]

1 1

10

100

Re

Oil Equation [6.3-3] Figure 6.16.

1000

10000

100000

Hydrogen Equation [6.3-4]

Correlations for the numerically obtained mass transfer data

6.3.7. Packed bed model The validation of the intra-particle concentration profiles and the verification of the kinetic and diffusional models proposed by Ramírez et al. (2004, 2006) allowed to create a fullyheterogeneous 3D packed bed reactor model in order to study the supercritical hydrogenation process in a semi-industrial scale. The packed bed model was created using a unit cell approach (a repetitive geometrical model) and simulations were carried by stages until full conversion was achieved for every case analyzed. Particle size of Dp = 0.47 mm was chosen for the simulations because of being considered the optimal particle size for the hydrogenation reaction, as explained in prior sections of this chapter. For every simulation stage, obtained contour fields for all computed variables (velocity and species concentration) were recorded at the model flow outlet, and fed as inlet boundary conditions for the following stage. This simulation strategy is highly timedemanding, but it guaranteed that the computational size of the CFD model was kept in reasonable sizes. Figure 6.17 shows an example of the obtained concentration profiles at different simulation stages along the packed bed reactor. For all cases analyzed, it can be noticed that there is a fast conversion in the first sector of the packed bed reactor. This fact is related to the extremely high reaction rates achieved in this zone of the reactor, approximately 1000 times higher than for the traditional biphasic hydrogenation process (Tacke et al., 1996). This fact can be clearly noticed when analyzing the conversion profiles along the packed bed length (see Figure 6.18). The high reaction rates obtained at the bed inlet are reflected in the small reactor sizes obtained (less than 10 cm length for Re = 200).

- 132 -

Chapter Six

Figure 6.17.

Linoleic fatty acid concentration contour fields for different stages of the packed bed reactor model (T= 473 K; Re = 100)

Iodine Value (IV)

150

100

50

0 0

0.02

0.04

0.06

0.08

0.1

Reactor length [m] Re = 10 Figure 6.18.

Re = 40

Re = 70

Re = 100

Re = 150

Re = 200

Conversion profile along the packed bed length for different inlet velocities at 457 K

- 133 -

Chapter Six

As it is important to minimize the presence of trans C18:1 (elaidic fatty acid) in the hydrogenation product due to its health implications (Ramírez, 2005), oil fatty acids mass fraction in the hydrogenated product was also obtained (see Figure 6.19). It can be seen that an increase in the reaction temperature makes the reaction move primarily to the formation of stearic fatty acid, reducing the formation of elaidic fatty acid. An increase in the reaction temperature produces an increase of the reaction rates, favoring the formation of stearic fatty acid.

Mass fraction

0.6

L

0.4

0.2

O S

E

0 120

100

E - 457 K E - 473 K Figure 6.19.

Iodine value (IV)

L - 457 K L - 473 K

80

O - 457 K O - 473 K

60

S - 457 K S - 473 K

Oil fatty acids mass fraction in the hydrogenated product (Re = 200)

CONCLUSIONS Computational Fluid Dynamics proves to be a useful tool when applied to the study of the catalytic hydrogenation of sunflower oil in the presence of a supercritical solvent. This computational technique allows to obtain in a fast and economical way information that otherwise is extremely complicated and expensive to obtain experimentally. Numerical simulations for the supercritical hydrogenation of sunflower oil were validated against experimental data previously published by this research group (Ramírez et al., 2006), obtaining a good agreement in every case analyzed. Magnitude of relative errors was evaluated for the developed simulations, finding that this value was within the magnitude order of the experimental error previously reported by our research group. Due to the difficulty for obtaining experimental data at very low particle sizes, a CFD simulation extrapolating the catalyst evaluation for Dp = 0.1 mm was also made in order to corroborate the catalyst effectiveness. It can be seen in Figure 6.9 that for a particle with Dp = 0.1 mm, numerically obtained concentration profiles were almost identical to those reported experimentally for Dp = 0.47 mm (Ramírez et al., 2006), reassuring the idea that for a particle diameter of 0.47 mm or lower, the effectiveness factor of the catalyst is close to the unity.

- 134 -

Chapter Six

Obtained numerical data allowed analyzing local mass transfer coefficients over the catalyst surface. Sherwood number for the process was also evaluated, and it was found that its value lies in between those expected for a gas and those expected for a liquid. Particle size, temperature and superficial velocity effects on external mass transfer were analyzed. It was found that mass transfer coefficients increased with an increase in temperature and in particle size. When analyzing the effects of the superficial velocity over the mass transfer coefficients, it was found that for high superficial velocities the system shows the behavior of a well-stirred reactor. It was observed that at high superficial velocities there is no effect over mass transfer coefficients. Numerically obtained Sherwood numbers were correlated as a function of Reynolds and Schmidt numbers, in order to present a simple set of equations capable to predict mass transfer coefficients in the supercritical hydrogenation of sunflower oil. Correlations presented are valid for predicting mass transfer in a range of 30 < Re < 12954 and 0.7 < Sc < 2.1, whit an AARD = 2.8 %. The validation of a single particle model allowed creating a 3D heterogeneous packed bed model in order to study the hydrogenation reaction in a semi-industrial scale, obtaining concentration and conversion profiles along the packed bed. Reaction rates obtained in the first 10 % of the length of the reactor are extremely high in every case analyzed (approximately 1000 times higher than for the traditional biphasic hydrogenation process). Further work must be done in the study of a packed bed reactor configuration in order to find optimal design and operation conditions for the industrial development of the supercritical hydrogenation process.

REFERENCES Bode, J., (1994). Applications of computational fluid dynamics in the chemical industry. Chemical Engineering & Technology, 17, 145 – 148. Brunner, G., (1985). Mass transfer in gas extraction, a supercritical fluid technology. In: J.M.L. Penninger, M. Radosz and M.A. McHugh (Eds.). Elsevier Science Publishers, Amsterdam. Dullien, F.A.L., (1979). Porous media: Fluid transport and pore structure. Academic Press, New York. Farrauto R.J., and Bartholomew C.H., (1997). Fundamentals of Industrial Catalytic Processes. London, UK: Chapman & Hall. Guardo, A., Casanovas, M., Magaña, I., Martínez, D., Ramírez, E., Larrayoz, M.A., Recasens, F., (2006). CFD modeling on external mass transfer and intra-particle diffusional effects on the supercritical hydrogenation of sunflower oil. Chemical Engineering Science (submitted). Harris, C.K., Roekaerts, D., Rosendal, F.J.J., Buitendijk, F.G.J., Daskopoulos, Ph., Vreenegoor, A.J.N., Wang, H., (1996). Computational fluid dynamics for chemical reactor engineering. Chemical Engineering Science, 51, 1569 – 1594. Härröd, M., Van den Hark, S., Macher, M.B., Moller, P., (2001). Hydrogenation at supercritical single-phase conditions. In: Bertucco, A., and Vetter, G., (Eds.). High Pressure Process Technology: Fundamentals and Applications. Elsevier Science Publishers, New York.

- 135 -

Chapter Six

Hashimoto, K., Muroyama, K., Nagata, S., (1971). Kinetics of the hydrogenation of fatty oils. Journal of the America Oil Chemists’ Society, 48, 291 – 295. Jackson, R., (1977). Transport in porous catalysts. Elsevier, New York. Kuipers, J. A. M., and Van Swaaij, W. P. M., (1998). Computational fluid dynamics applied to chemical reaction engineering. In: Wei, J., Anderson, J., Bischoff, K., Denn, M., Seinfield, J., Stephanopoulos, G., (Eds.); Advances in Chemical Engineering 24; Academic Press: San Diego, CA; pp. 227 – 328. Petersen, E.E., (1965). Chemical reaction analysis. Prentice-Hall, New Jersey. Poling, B., Prausnitz. J., O’Connell, J., (2000). The properties of gases and liquids. McGraw-Hill, New York. pp. 3.1 – 10.56. Ramírez, E., Recasens, F., Fernández, M., Larrayoz, M.A., (2004). Sunflower oil hydrogenation in Pd/C in SC propane in a continuous recycle reactor. AIChE Journal, 50, 1545 – 1555. Ramírez, E., (2005). Contribution to the study of heterogeneous catalytic reactions in SCFs: Hydrogenation of sunflower oil in Pd catalysts at single-phase conditions. Ph.D. Thesis. Chemical Engineering Department, Universitat Politècnica de Catalunya, Barcelona. Ramírez, E., Larrayoz, M.A., Recasens, F., (2006). Intraparticle diffusion mechanisms in SC sunflower oil hydrogenation on Pd. AIChE Journal, 52, 1539 – 1553. Ranade, V., (2002). Computational flow modeling for chemical reactor engineering. Academic press, New York, pp. 244 – 422. Ranz, W.E., and Marshall Jr., W.R., (1952). Evaporation from drops, part 1. Chemical Engineering Progress, 48, 173 – 180. Reid, R., Prausnitz, J., Poling, B., (1987). The properties of gases and liquids. McGraw-Hill, Boston. pp. 95 – 205. Satterfield, C.N., (1970). Mass transfer in heterogeneous catalysis, 2nd ed., McGraw-Hill, New York. Tan, C.-S., Liang, S.-K., Liou, D.-C., (1988). Fluid-solid mass transfer in a supercritical fluid extractor. Chemical Engineering Journal, 38, 17 – 22. Tacke, T., Wieland, S., Panster, P., (1996). Hardening of fats and oils. In: Rudolf von Rohr, P., and Trepp, Ch., (Eds.); Proceedings of the 3rd International Symposium on High-Pressure Chemical Engineering. Zurich, Switzerland. Elsevier Science Publishers; pp. 17 - 22. Tacke, T., Wieland, S., Panster, P., (2003). Selective and complete hydrogenation of vegetable oils and free fatty acids in supercritical fluids. In: De Simone, J.M., (Ed.) Green Chemistry Using Liquid and Supercritical Carbon Dioxide. Oxford, UK: Oxford Univ. Press. Wheeler, A., (1951). Advances in catalysis, Vol. 3. Academic Press, New York. Wheeler, A., (1955). In: P.H. Emmett (Ed.), Catalysis, Vol. 2. Reinhold, New York. White, F., (1991). Viscous fluid flow. McGraw-Hill, New York. pp. 482 – 542.

- 136 -

Chapter Six

Yaws, C., (1999). Chemical properties handbook. McGraw-Hill, New York. pp. 1 – 55. Youngquist, G.R., (1970). SYMPOSIUM ON FLOW THROUGH POROUS MEDIA: Diffusion and flow of gases in porous solids. Industrial & Engineering Chemistry, 62 (8), 52 – 63.

- 137 -

CHAPTER SEVEN

conclusions and future work

Conclusions

CFD proves to be a useful tool when modeling mass and heat transfer and chemical reactions in packed bed equipment. It allows obtaining velocity, pressure, temperature and species concentration contour plots inside the bed, reaching information almost impossible to get by means of experimental work. Using CFD commercially available codes (a Finite Element code, Comsol 3.x, and a Finite Volume code, Fluent 5.x/6.x) allowed us to simulate and validate velocity and temperature profiles around simple geometries, and to study mass and heat transfer phenomena and heterogeneous reaction at low and high pressure in packed beds. A summary of the conclusions reached in each section of this thesis work is presented in this chapter.

7.1.

Validation models

Two computational flow models were developed to validate flow and heat transfer around spheres (single or stacked). For the flow validation test, flow around a simple cubic stack of spheres was used as validation model to test the capabilities of the solver reproducing experimental data. The model predictions were verified by comparing the simulation results with the published experimental and computational results. Predicted results showed excellent agreement with the experimental data of Suekane et al., (2003). Mesh sensitivity was established and optimal average mesh density for flow problems could be obtained. For the heat transfer validation test, a sphere suspended in an infinite fluid was used as validation tool to test the capabilities of the solver reproducing an analytical solution. Drag coefficient over particle surface was recorded and compared against the Stokes’ law and the graphical correlation presented by Lapple and Shepherd (1940), obtaining an overall good agreement between the compared sets of data, and a neglectable mesh dependency on the results. In the case of the prediction of heat transfer parameters, mesh sensitivity tests were performed, and optimal average mesh density over the heat transfer surface was established. Numerical results obtained were compared against the theoretical answer for estimating the heat transfer coefficient obtained by Ranz and Marshall (1952), obtaining a good agreement between numerical and theoretical answers.

7.2.

Wall-to-fluid heat transfer

CFD proves to be useful in the wall to fluid heat transfer parameter estimation, and also for calculation of pressure drop along the bed in packed bed equipment. It was possible to model a realistic case of a packed bed of spheres including contact points within the surfaces involved in the geometry. The calculated velocity profiles fitted qualitatively the expected results, and the calculated values of pressure drop along the bed adjust quite well with previously published and accepted correlations. Flow structures within the bed (i.e., wall channeling, stagnant points, eddy flows) were easily identifiable. Obtained temperature profiles inside the bed allowed estimating wall heat transfer parameters such as Nuw and kr/kf . The results obtained for all of the cases studied (laminar and turbulent) agree among themselves and with the selected empirical and semi-empirical correlations when analyzing pressure drop and effective radial conductivity. This can be explained by the similarity in the

- 141 -

Conclusions

velocity field obtained for each simulation. The calculation of these parameters is more closely related to velocity fields than to mixing parameters. The prediction of the mixing rate within the bed along with the near-wall treatment appreciably affects the estimate of Nuw. Flow regime zones can be identified using the heat transfer coefficient estimate. A laminar solution overestimates the value of this coefficient in the turbulent flow zone, and turbulent solutions tend to underestimate the value of the coefficient in the laminar transition zone. Turbulence models used for simulations do not predict the transition regime well, and this can be seen in the discrepancy between numerically obtained results and correlations in the low Re and transition range. The definition of a good mesh allows calculations of fluid dynamics variables, as velocity and pressure. However, in our case, the proposed geometry governs mesh density and element size in the near-wall area, and this fact affects the definition of an appropriate y+ parameter in order to apply a correct near-wall treatment for certain turbulence models (e.g. the κ − ε family models). To define an adequate y+ for a correct coupling for the two-equation models at the near-wall treatment is not an easy task. Therefore, the y+ parameter is crucial during the selection of the appropriate turbulence model to apply within the simulation. A good near-wall modeling is fundamental to obtain more accurate results in pressure drop and heat transfer calculations and the selection of the right turbulent model will depend on the geometry proposed and the values obtained for y+ in the wall. Results obtained with the Spalart–Allmaras turbulence model show better agreement than the two-equation RANS models for pressure drop and heat transfer parameter estimation. This could be due to the fact that this model uses a coupling between wall functions and damping functions for near-wall treatment, does not include additional diffusion or dissipation terms in its formulation and does not present the stagnation point anomaly. Factors such as the misestimating of ε, κ or µt can lead to differences in flow and temperature profiles that can be translated into miscalculation of heat transfer parameters. Therefore, the Spalart – Allmaras model could be a good tool for these kinds of flows because the y+ problem is solved automatically in spite of the necessity of checking the upper values of y+ (less than 120). Turbulence models used for simulations do not properly predict the transition regime and this can be observed in the discrepancy between numerically obtained results and correlations in the low Re and transition range. Results present slower ratios of convergence at low Re than at high Re.

7.3.

Particle-to-fluid convective mass and heat transfer at low and high pressure

CFD proves to be a reliable tool when modeling convective mass and heat transfer phenomena in packed beds. It allows to analyze either free or forced convection situations and the obtained results can be compared qualitatively and quantitatively against previous published data. Forced convection at low pressure in a packed bed was simulated, and the influence of velocity over mass and heat transfer could be analyzed. For values of Re > 10 it was obtained a good agreement with the correlations presented by Wakao and Kaguei (1982) for mass transfer, and Wakao et al. (1979) for heat transfer. At lower Reynolds numbers (Re < 10), results shows that the fitting against correlation is not good. This situation is especially notorious when analyzing heat transfer. For a single velocity condition different meshes give results in a wide range of Nu and no relation with mesh density can be established in any case. No mesh sensitivity was

- 142 -

Conclusions

noticed for the laminar and transition flow zones, but in the turbulent flow zone a good definition of the mesh around the particles surface is of primal importance in order to capture the turbulence vorticity energetic scales associated effects. Mixed (free+forced) convection at high pressure in a packed bed was also analyzed. For a supercritical fluid in laminar flow regime it was possible to study the effects of the density gradient, flow direction and velocity over mass and heat transfer. It was noticed that the presence of large density gradients conditioned the mesh influence over the numerical results when computing the mixing and the mass/heat transfer within the computed domain. An adverse density gradient generates hydrodynamical instabilities that produce an increase of the axial dispersion and a diminished mass/heat transfer rates. This fact caused numerical instability in the simulation process. In order to eliminate numerical instability, the mesh was locally refined in each case trying to avoid the presence of high density and velocity gradients along the bed. Influence of flow direction over mass/heat transfer was also analyzed and it was noticed that in assisting flow regime, greater heat transfer rates are obtained. The rising in the density along the bed height due to the adverse density gradient helps the axial dispersion to grow, obtaining less mixing within the bed and smaller heat transfer coefficients. Flow velocity also affects the heat transfer rate. The value of Sh and Nu increases almost linearly with flow velocity and this increase is due exclusively to the increase in the contribution of forced convection to the overall mass/heat transfer. Free convection is independent of flow velocity and its value remains almost constant within the studied velocity range. A novel correlation (Eq. [5.3-3]) is presented for estimating the free and forced convection effects and Nu from Re, Pr and GrH . Presented correlations are valid for 9 < Re < 96, 2.2 < Pr < 3.3, and 1 × 108 < GrH · Pr < 4 × 1010. The obtained numerical results validates the idea that the modified correlation presented (Guardo et al., 2006), based on the one proposed by Stüber et al. (1996) for mixed convection mass transfer under supercritical conditions in packed beds, can be used to describe heat transfer phenomena in a supercritical packed bed reactor under mixed convection regime at high pressures.

7.4.

High pressure heterogeneous reaction

Computational Fluid Dynamics proves to be a useful tool when applied to the study of the catalytic hydrogenation of sunflower oil in the presence of a supercritical solvent. This computational technique allows to obtain in a fast and economical way information that otherwise is extremely complicated and expensive to obtain experimentally. Numerical simulations for the supercritical hydrogenation of sunflower oil were validated against experimental data previously published by this research group (Ramírez et al., 2006), obtaining a good agreement in every case analyzed. Magnitude of relative errors was evaluated for the developed simulations, finding that this value was within the magnitude order of the experimental error previously reported by our research group. Due to the difficulty for obtaining experimental data at very low particle sizes, a CFD simulation extrapolating the catalyst evaluation for Dp = 0.1 mm was also made in order to corroborate the catalyst effectiveness. It can be seen in Figure 6.9 that for a particle with Dp = 0.1 mm, numerically obtained concentration profiles were almost identical to those reported

- 143 -

Conclusions

experimentally for Dp = 0.47 mm (Ramírez et al., 2006), reassuring the idea that for a particle diameter of 0.47 mm or lower, the effectiveness factor of the catalyst is close to the unity. Obtained numerical data allowed analyzing local mass transfer coefficients over the catalyst surface. Sherwood number for the process was also evaluated, and it was found that its value lies in between those expected for a gas and those expected for a liquid. Particle size, temperature and superficial velocity effects on external mass transfer were analyzed. It was found that mass transfer coefficients increased with an increase in temperature and in particle size. When analyzing the effects of the superficial velocity over the mass transfer coefficients, it was found that for high superficial velocities the system shows the behavior of a well-stirred reactor. It was observed that at high superficial velocities there is no effect over mass transfer coefficients. Numerically obtained Sherwood numbers were correlated as a function of Reynolds and Schmidt numbers, in order to present a simple set of equations capable to predict mass transfer coefficients in the supercritical hydrogenation of sunflower oil. Correlations presented are valid for predicting mass transfer in a range of 30 < Re < 12954 and 0.7 < Sc < 2.1, whit an AARD = 2.8 %. The validation of a single particle model allowed creating a 3D heterogeneous packed bed model in order to study the hydrogenation reaction in a semi-industrial scale, obtaining concentration and conversion profiles along the packed bed. Reaction rates obtained in the first 10 % of the length of the reactor are extremely high in every case analyzed (approximately 1000 times higher than for the traditional biphasic hydrogenation process). Further work must be done in the study of a packed bed reactor configuration in order to find optimal design and operation conditions for the industrial development of the supercritical hydrogenation process.

7.5.

Future work

Further work on the study of heat and mass transfer phenomena in packed beds can be done applying a CFD methodology, especially when dealing with supercritical fluids. The advantages of CFD simulations (low cost, fast and reliable information) can be applied to situations where obtaining experimental data is technically complex and not cost effective. When dealing with supercritical fluids, materials and equipment needed to deal with the high pressures and high temperatures required for the experimentation are extremely expensive, limiting sometimes the experimentation possibilities. Further work can be done in this subject, using CFD to simulate complex situations involving supercritical fluids, as radial effects in packed beds or axial dispersion studies when natural convection is involved into the heat/mass transfer phenomena. Further research can be done for studying the effects of the configuration of the packing structure in the efficiency of the transport process, and the effects of the geometrical configuration (particle size, pore size and pore distribution) on transport and reaction parameters.

REFERENCES Guardo, A., Coussirat, M., Recasens, F., Larrayoz, M.A., Escaler, X., (2006). CFD study on particle-to-fluid heat transfer in fixed bed reactors: Convective heat transfer at low and high pressure. Chemical Engineering Science, 61, 4341 – 4353.

- 144 -

Conclusions

Lapple, C.E., and Shepherd, C.B., (1940). Calculation of particle trajectories. Industrial & Engineering Chemistry, 32, 605 – 617. Ramírez, E., Larrayoz, M.A., Recasens, F., (2006). Intraparticle diffusion mechanisms in SC sunflower oil hydrogenation on Pd. AIChE Journal, 52, 1539 – 1553. Ranz, W.E., and Marshall Jr., W.R., (1952). Evaporation from drops, part 1. Chemical Engineering Progress, 48, 173 – 180. Suekane, T., Yokouchi, Y., Hirai, S., (2003). Inertial flow structures in a simple-packed bed of spheres. AIChE Journal, 49, 10 – 17. Stüber, F., Vázquez, A.M., Larrayoz, M.A., Recasens, F., (1996). Supercritical fluid extraction of packed beds: external mass transfer in upflow and downflow operation. Industrial & Engineering Chemistry Research, 35, 3618–3628. Wakao, N., Kaguei, S., Funazkri, T., (1979). Effect of fluid dispersion coefficients on particle-tofluid heat transfer coefficients in packed beds: correlation of Nusselt numbers. Chemical Engineering Science, 34, 325 – 336. Wakao, N., and Kaguei, S., (1982). Heat and mass transfer in packed beds. Gordon and Breach Science Publishers, New York. pp. 138 – 160; 264 – 295.

- 145 -

APPENDICES

APPENDIX A

BIBLIOGRAPHIC REVIEW

on CFD design issues, supercritical fluids extraction and reaction

Luke Andrews, « Open book »

Michael Graf, « Turbulence »

Appendix A

A-1. Extraction of natural products A wide range of extractions of natural products from vegetal, animal and wood/fibre material can be carried out with a supercritical solvent (see Tables A-1 to A-3), and there is a growing interest in developing such processes, specially for the food, cosmetics and pharmaceutical industries. Extract

T [K]

P [M Pa]

Solvent

Oil

313

12

CO2

Australian ginger root Black pepper (Piper nigrun L.)

Oil Essential oil

282 – 308 303 – 323

6 – 10 15 – 30

Borage (Borago officinalis L.) seeds

Total extract

313

10 – 35

CO2 + ethanol CO2 CO2 + caprylic acid methyl ester

305 – 348

7.5 – 30

CO2

308

20 – 30

Oil

315

30

Total extract

323 – 343

25 – 35

CO2 CO2 + ethyl alcohol CO2 Ethane

308 – 338

15 – 45

CO2

Simándi et al., 2002

293 – 323

8 – 28

CO2

Mira et al., 1999

Total extract

313 – 323 295 – 323 311 – 323 313 – 333

25 – 30 7 – 30 20 – 30 20 – 30

CO2 CO2 + ethanol CO2 + ethanol CO2 + ethanol

Catchpole at al., 2002

Total extract

303 – 313

10 – 20

CO2

Povh et al., 2001

β - carotene

330

25

Subra et al., 1998

Oleoresin

303 – 313

20 – 25

CO2 CO2 + ethanol CO2 + propanol

333 – 393

24 – 31

CO2

Chiu et al., 2002

293 – 313 313 – 343

7 – 25 10 – 40

CO2 CO2 + H2O CO2 Propane

Kiriamiti et al., 2003 Saldaña et al., 2002

CO2

del Valle et al., 2003

Raw material Angelica archangelica L. root

Caraway (Carum carvi L.) seeds Coriander seed Corn germ Cupuaçu (Theobroma grandiflorum).,Brazilian fruit Dandelion leaves Dehydrated orange peel Dried Saw Palmetto berries St. John’s Wort flowers Echinacea purpurea (aerial parts) Kava roots and stems Flowers of chamomile (Chamomilla recutita [L.] Rauschert) Freeze dried carrots Ginger (Zingiber officinale Roscoe)

Carvone, Limonene Oil

β - Amyrin β - Sitosterol Oil Limonene

Ground pyrethrum flowers Guaraná seeds

Terpene lactones and flavonoids Pyrethrins Caffeine

Hiprose fruit

Total extract

298 – 308

5 – 25

Oleoresin

313

12 – 32

Ginkgo leaves

Jalapeño peppers (Capiscum annuum L.) prepelletised Lemon balm (Melissa Officinalis L.) Lemongrass (Cymbopogon citratus) Lovage (Levisticum officinale Koch.) Nutmeg (Myristica fragans Houttuyn) Pungent spice paprika powder Shells of the bacuri fruit (Platonia insignis Mart) Soybeans Summer savory (Satureja hortensis L.)

Table A-1.

Total extract Antioxidants Essential oil Essential oil

References Doneanu and Anitescu, 1998 Badalyan et al., 1998 Ferreira et al., 1999 Daukšas et al., 2002 Baysal and Starmans, 1999 Illés et al., 2000 Rónyai et al., 1998 de Azevedo et al., 2003

Zancan et al., 2002

Illés et al., 1997

308 – 313

10 – 18

CO2

Ribeiro et al., 2001

296 – 323 313 – 323

8.5 – 12 8 – 35

CO2 CO2

Carlson et al., 2001 Daukšas et al., 1999

Oil

296

9

CO2

Spricigo et al., 1999

Oleoresin

308 – 328

10 – 40

CO2 Propane

Daood et al., 2002

Total extract

289 – 323

6.3 – 20

CO2 + ethanol

Monteiro et al., 1997

Oil

333 – 353

16 – 69

CO2 + ethanol

Montanari et al., 1999

Oil

313

12 – 18

CO2

Esquível et al., 1999

Vegetal products extractions carried out in SCF

- 151 -

Appendix A

Extract

T [K]

P [M Pa]

Solvent

References

Oil

308

34.5

CO2

Dunford et al., 1998

Black scabbardfish (Aphanopus carbo)

PCB, chlorinated pesticides

309 – 337

10 - 24

CO2

Antunes et al., 2003

Cow brain

Cholesterol

323 - 343

23 – 27

CO2

Dried egg yolk

Phospolipids

313

51.7

CO2

Pigskin Raw wool

Fat Natural wax Fat

10 - 60 7 - 20 10 – 24 10.4 – 24

CO2 CO2 + ethanol

Sheepskin

293 – 313 303 - 353 318 313 353

25

CO2 + toluene

Raw material Atlantic mackerel (Scomber scombrus)

Wool scour wastes

Table A-2.

Wax

Fish and animal material extractions carried out in SCF

Raw material

Extract

Cellulose-based filter papers

Sugar cane bagasse Pinus Taeda wood chips Wood pulp

Table A-3.

CO2

Vedaraman et al., 2005 Boselli and Caboni, 2000 Vaquero et al., 2006 Eychenne et al., 2001 Marsal et al., 2000a Marsal et al, 2000b López-Mesas et al., 2005

T [K]

P [M Pa]

Uranyl nitrate

333

25

Lignine

415 – 471

14.7 - 23

Total extract

353

51.6

Solvent CO2 + methanol, acetyl acetone, methylisobutyl ketone, Tributyl phosphate CO2 + ethanol/water CO2 + methanol

References

Shamsipur et al., 2001

Pasquini et al., 2005 McDaniel et al., 2001

Wood and fiber material extractions carried out in SCF

A-2. Soil remediation Decontamination of soils using SCF is an attractive process compared to extraction with liquid solvents because no toxic residue is left in the remediated soil and in contrast to thermal desorption soils are not burned. Specially, the removal of typical industrial wastes (PAH, PCB, fuels) can be easily achieved through SCE (see T able A-4). Contaminant PCB 1,1,1-trichloroethane (DDT) Toxaphene

Removal efficiency

T [K]

P [M Pa]

Solvent

References

70 – 90 %

313

10

CO2

Brady et al., 1987

Polycyclic aromatic hydrocarbons (PAH)

~ 100 %

313 – 423

40

Carbamate pesticides

39.6 – 91.7 %

327

37.8

> 88 %

353

34

Gasoline Diesel

Table A-4.

CO2 CO2 + methanol Acetone Hexane Methanol Toluene CO2 CO2 + methanol CO2

Examples of applications of SCE to soil remediation

- 152 -

Barnabas et al., 1995 Lojková et al., 2005

Izquierdo et al., 1996 Yang et al., 1995

Appendix A

A-3. Heterogeneous catalysis The use of SCF as a reaction media can be a real advantage when using heterogeneous catalysts, since the diffusion rates are enhanced compared to reactions in the liquid phase. See Table A-5 for further details on chemical reactions carried out in SCF. T [K] ALKYLATION

P [MPa]

References

323 – 473

3.5 – 5

Fan et al., 1997a

CO2

323 – 413

3.4 – 15.5

Clark and Subramaniam, 1998

Deloxan

Propene CO2

433 – 573

15 – 20

Hitzler et al., 1998

Zeolite

Heptane

3.4

Dardas et al., 1996

593 – 598

3.36 – 5.6

Collins et al., 1998

Catalyst type

Solvent

Isopentane and isobutene Isobutane and isobutene

Zeolite

Isopentane Isobutane

1-butene and isobutane

Zeolite

Reaction

Mesytilene and propene Mesytilene and propan-2-ol

CRACKING Heptane

598

DISPROPORTIONATION Toluene to p-xylene and benzene 1,4-diisopropylbenzene to cumene and 1,3,5-triisopropylbenzene Ethylbenzene to benzene and diethylbenzene

Zeolite

Benzene n-pentane

533

20

Zeolite

Butane Pentane

573 – 673

5

Oleic acid and methanol

Sulfonic

Zeolite

Tiltscher et al., 1984 Tiltscher and Hofmann, 1986 Niu and Hoffmann, 1995, 1997

ESTERIFICATION CO2

> 823

0.95 – 1.3

Vieville et al., 1993

FISHER – TROPSCH SYNTHESIS Ru/Al2O3 Co/SiO2

n-hexane

313 – 341

4.5

Fe

n-hexane

513

7–8

1-hexene

Pt/Al2O3

CO2+n-pentane CO2+n-hexane

1-hexene to olefinic oligomers

Pt/Al2O3

CO and H2 to liquid hydrocarbons

Yokota et al., 1990, 1991a, 1991b Snavely and Subramaniam, 1997

ISOMERIZATION 291

18

300.7

27.7

Clark and Subramaniam, 1996 McCoy and Subramaniam, 1995

OXIDATION Toluene to benzaldehyde Isobutane to tert-butyl alcohol

Co/Al2O3

CO2

Pd/C

281

8

Dooley and Knopf, 1987

426

4.4 – 5.4

Fan et al., 1997b

HYDROGENATION Fats and oils Acetophenone Cyclohexene Cyclohexene Unsaturated ketones

Deloxan

CO2 - Propane

323 – 433

8 – 16

Tacke et al., 2003

Deloxan

CO2 - Propane

313 – 593

6 – 12

Hitzler and Poliakoff, 1997

Pd/C Pd/Al2O3

CO2 CO2

343 323 – 493

13.6 12 – 17.5

Cu

Propane

473 – 573

15

α - pirene

Pd/C

CO2

323

14

Palm oil

Pd/C

Propane

338 - 408

15

Ni

CO2 Propane Dimethyleter

393 - 413

14

428 - 488

20

Arunajatesan et al., 2001 Bertucco et al., 1997 van den Hark et al., 1999, 2000, 2001 Chouchi et al., 2001 Macher and Holmqvist, 2001 King et al., 2001 Ramirez et al., 2004 Ramirez, 2005

Fatty acid methyl esters

Soybean oil Sunflower oil

Table A-5.

Pd/C

Survey of heterogeneous catalytic reactions carried out under supercritical conditions (Baiker, 1999; Ramirez et al., 2002)

- 153 -

Appendix A

A-4.

CFD in chemical engineering process equipment design

and

Typical process equipment, key design issues, and perspectives on the current status of CFD for simulating flows of practical interest in chemical engineering are summarized in Tables A-6 to A-11. (Joshi and Ranade, 2003). Typical equipment Design issues / applications

Pipes, pipe fittings, T mixers, valves, flow meters, membrane modules, heat exchangers, monoliths, etc. Pressure drop, heat transfer coefficient, mixing, residence time distribution, identification of hot spots/dead regions, influence of non-uniform material properties, extrapolation to non-standard configurations, end effects.

CURRENT STATUS Knowledge of physics

CFD codes / design applicability

Excellent for laminar flow of simple fluids. Constitutive equations for complex fluids are not well developed. Reasonable for turbulent flow of simple fluids. Turbulence models for non-Newtonian fluids are not well developed. Commercial CFD codes/models allow implementation of different constitutive equations/non-uniform material properties. Grid quality, constitutive equations and turbulence models control the quality of simulations. It is essential to ensure that the appropriate turbulence model is selected for simulating complex turbulent flows even for qualitative screening of configurations; reasonable guidelines for such a selection are available. Applicability is generally excellent for laminar flows of Newtonian and simple nonNewtonian fluids; limited for viscoelastic/structured fluids. CFD simulations are widely used for addressing key design issues.

PATH FORWARD Limitation of current models

Geometrical complexity: grid quality, curvature effects. Rheological complexity: constitutive equations, solvers/algorithms. Turbulence: transition regime, transport rates near walls, computing resources. Reactive mixing: computing resources, closure models.

Experiments

Physical models

Influence of small-scale features performance (Liou et al., 2002).

Development needed for overcoming the limitations

on

Behavior of complex fluids in nonrheometric flows (Rothstein and McKinley, 2001). Data at the onset of turbulence (Palikaras et al., 2002). Measurement of wall transport rates for model discrimination (Vogel, 1984).

Influence of curvature and severe changes in flow direction on turbulence or complex fluids (Suga, 2003). Constitutive equations based on molecular structure (van Ruymbeke et al., 2002). Turbulence models adaptable to the resolution requirements (Speziale, 1998). Understanding the role of unsteady flows (spatially or temporally) in mixing and other transport processes (Jana et al., 1994).

Measurement of local concentration and segregation profiles (Verschuren et al., 2002).

Table A-6.

Typical process equipment, key design issues, and perspectives on the current status of CFD for simulating single-phase flows through conduits/channels (Joshi and Ranade, 2003)

- 154 -

Appendix A

Typical equipment Design issues / applications

Knowledge of physics

CFD codes / design applicability

Limitation of current models

Development needed for overcoming the limitations

Table A-7.

Falling film evaporators, fiber spinning, spray coating, multiphase flow through monoliths, etc. Shape of free surface/film thickness, influence of surface characteristics (roughness/wetting), non-uniform material properties, heat/mass transfer coefficients, hot spots, learning tool for developing macroscopic closure models. CURRENT STATUS Basic models to capture the role of surface forces on the shape of the interface between immiscible phases are available. Understanding of surface characteristics, contact angle, and wetting and wall adhesion is not adequate. Quantitative prediction of instabilities and regime transitions is not yet adequately possible. Commercial CFD codes allow simulating free surface flows with the “volume of fluid” (VOF) approach or other front tracking methods. Demands on grid resolution/computing resources are high: usually threedimensional simulations are needed; small time step and several iterations per time step are needed for simulating transients accurately; convergence behavior is sensitive to the grid quality VOF-based codes are relatively robust compared to those with surface tracking; simulations require expert intervention. CFD simulations are used more as a learning tool than as a design tool; in some instances, calibrated models may be used for design. PATH FORWARD Wetting and wall adhesion: hysteresis in contact angle, grid quality, curvature effects. Transport at interface: capturing small-scale features of the interface, computation of gradients at the interface. Flow regimes: grid quality, transport of the interface, computational resources. Prediction of “die-swell” and other peculiar characteristics of viscoleastic fluids. Experiments Physical models Characterization of wetting and drop Models for wetting characteristics dynamics on different surfaces (Ted (Fukai et al., 1995) Mao et al., 1997). Influence of gradients of surface Data on mass and heat transfer through tension on transport rates the interface (Gulawani, 2002). Characterization of film flows and Capturing small-scale structures on associated instabilities (Ambrosini et al., the interface. 2002). Data on the influence of viscoelasticity Interface models to capture regime on evolution of free surfaces (Coopertransitions. White et al., 2002). Constitutive equations suitable for practical flows of interest. Typical process equipment, key design issues, and perspectives on the current status of CFD for simulating free surface flows (Joshi and Ranade, 2003)

- 155 -

Appendix A

Typical equipment Design issues / applications

Rotating disk contactors, stirred reactors/mixers, pumps/turbines, screw extruders, rotary kilns/dryers, ball mills, centrifuges, etc. Number, shape, size, location and speed of impellers or rotating walls, mixing/residence time distribution, heat and mass transfer, scale-up/scale-down

CURRENT STATUS Adequate for single-phase flows in stirred vessels.

Knowledge of physics

Models for simulating dispersion or suspension of multiphase flows are inadequate. Models for granular flows in kilns/ball mills are inadequate. Commercial CFD codes/models allow simulations of rotating elements with either quasi-steady or sliding mesh approaches. Fine resolution is needed to capture vortices around rotating elements; full transient simulations based on the sliding mesh approach are very computation intensive.

CFD codes / design applicability

Prediction of mean flow is reasonably good. Quantitative prediction of turbulence characteristics is still not good. A priori prediction of multiphase flows in stirred vessels or rotary kiln is not possible; regime-specific models with calibrated parameters may provide useful design information. CFD simulations are used for understanding different design configurations mainly to assist engineering decision making.

PATH FORWARD Limitation of current models

Free surface at the top/vortex/gas entrainment. Turbulence: transition regime, non-Newtonian fluids. Low-frequency unsteady flows. Multiphase flows: coalescence/agglomeration, break-up, accumulation behind rotating elements, heat and mass transfer.

Experiments

Physical models Calibration of parameters of the turbulence models.

Development needed for overcoming the limitations

Data near the top free surface (Desai and Joshi, 1996). Measurement of turbulence kinetic energy, stresses and energy dissipation rate (Joshi et al., 2001) Characterization of unsteady flow and its impact on the performance (Kresta and Wood, 1993). Characterization of cavities formed behind rotating blades (Lu and Ju, 1989). Local measurements of phase fractions, phase velocities, and particle sizes at higher phase fractions (Bombač et al., 1997). Local measurements of heat transfer coefficients at different coil/jacket configurations (Karcz, 1999). Measurements of angle of internal friction, coefficient of restitution, granular pressure, and temperature for solidliquid flows.

Influence of the free surface/vortex on flow in the stirred vessels (Ciofalo et al., 1996) Turbulence models/wall functions for nonNewtonian fluids (Son and Singh, 2002). Ways of capturing low-frequency flow structures. Modified wall functions to capture heat transfer near walls with complex curvatures (Churchill and Zajic, 2002). Interface closures and coalescence/breakup models for high phase factions. Improved models for estimating the frictional viscosity of solid particles sliding on each other. Models for simulating formation of gas cavities behind the blades

Table A-8.

Typical process equipment, key design issues, and perspectives on the current status of CFD for simulating stirred or rotating vessels (Joshi and Ranade, 2003)

- 156 -

Appendix A

Typical equipment

Bubble/slurry columns fluidized beds, spray columns/extraction columns, cyclones, settling tanks, crystallizers/precipitators, distillation trays, condensers/evaporators, furnace/boilers, mist eliminators.

Design issues / applications

Distributor design, flow regime, minimal fluidization/suspension velocity, sensitivity with hardware changes, distribution of dispersed phase and interfacial area, mass and heat transfer, entrainment, scale-up, erosion.

CURRENT STATUS Reasonable for dilute multiphase flows.

Knowledge of physics

Inadequate for dense multiphase flows; closures for interface coupling for turbulent, dense dispersions are not adequate. Models involve several parameters, which have to be obtained from experimental data. Eulerian-Lagrangian and Eulerian-Eulerian approaches are used for simulating diluted and dense flows respectively; several particle trajectories need to be simulated for an adequate picture.

CFD codes / design applicability

A priori prediction of flow regimes is still not possible. Models are often calibrated without obtaining grid-independent results; ad hoc models to simulate the influence of turbulence or other particles are used. Though basic solvers for Eulerian-Lagrangian and Eulerian-Eulerian approaches are available in commercial CFD codes, user-defined functions are needed to implement problem-specific sub-models and closures. Use in practical design is often indirect: for qualitative screening and for assisting engineering decision making rather than for quantitative simulations.

PATH FORWARD Limitation of current models

Trajectory simulations in turbulent flows: drag coefficient, lift coefficient and virtual mass coefficient, turbulent dispersion of particles. Turbulence caused by dispersed particles. Momentum transport between different dispersed phases. Agglomeration/coalescence/break-up process. Phase change (evaporation/condensation).

Experiments

Physical models

Unambiguous experimental data to determine drag, lift and virtual mass coefficients under different conditions (Tomiyama et al., 2002).

Development needed for overcoming the limitations

Table A-9.

Experiments allowing separation of dispersion effects from lift and other effects. Data on Lagrangian quantities (Cassanello et al., 2001). Experiments to quantify the fraction of energy dissipated in the bulk of the continuous phase (away from the interface). Local measurements of size distribution to track evolution of it within the vessel (Buwa and Ranade, 2002). Measurement of the interfacial area and local mass transfer coefficients (Sun et al., 2003).

Direct numerical simulation (lattice Boltzmann or surface tracking) of multiphase flows to provide theoretical basis to empirical correlations of drag, lift and virtual mass transfer coefficients (Sankaranarayanan and Sundaresan, 2002). Models to ensure correct Lagrangian correlations (Gouesbet and Berlemont, 1999). Improved models for the interface momentum, energy and mass transfer. Better understanding of the formation of agglomeration/clusters and their influence of transport rates (Harris et al., 2002). Models to estimate the effective dispersion coefficients as a function of particle size and local turbulence length scales. Guidelines to select appropriate coalescence/break-up kernels.

Typical process equipment, key design issues, and perspectives on the current status of CFD for simulating multiphase flows (Joshi and Ranade, 2003)

- 157 -

Appendix A

Typical equipment Design issues / applications

Reactors for fast reactions (comparable time scales of mixing and reactions), combustors, burners, etc. Selectivity and yield, complete combustion and energy efficiency, NOx formation, nozzle designs, etc.

CURRENT STATUS Reasonable for single phase laminar flow.

Knowledge of physics

Closures for turbulent reactive mixing are available only for simple reaction schemes. Inadequate models for interfacial area limit modeling of multiphase reactive flows. Single phase reactive mixing models may be solved using the existing framework with additional closure models (probability density function/phenomenological approaches); user defined functions are needed for implementing specific phenomenological closures.

CFD codes / design applicability

Particle level reactions may be simulated with the Eulerian-Lagrangian approach. Simulations of reactive flows are computationally very demanding, and it is difficult to eliminate the influence of numerical issues completely. CFD simulations are used to address design issues of equipment involving single phase or dilute multiphase reactive flows. Capability of simulating multiphase reactive systems with the Eulerian-Eulerian or VOF approach is still in the primitive stages.

PATH FORWARD Limitation of current models

Development needed for overcoming the limitations

Table A-10.

Laminar reactive mixing in complex fluids. Turbulent reactive mixing in liquids. Dilute/dense multiphase reactive flows.

Experiments

Physical models

Measurements of reduction in striation thickness with time under wellcharacterized flows for complex fluids (Fourcade et al., 2001).

Better understanding of interaction between turbulent, mixing and reaction for cases with one of the reactants or products existing in a different immiscible phase (Marchisio et al., 2001).

Experiments to quantitatively measure the reduction in scale of mixing and intensity of mixing (Buchmann and Mewes, 2000). Measurements of cross correlations of concentrations for improved closures. Measurements of local interfacial area and transport coefficients from agglomerated particle clusters (van der Bosch, 1998). Measurements of transport coefficients under evaporation or condensation. Single particle combustion data under controlled flow situations.

Quantification of the role of complex rheology on stretching and folding to create large interfacial areas. Closure models for simulating complex reaction schemes. Models/correlations for interface transport rates from drops/particles with phase change. Models of cluster formation and to predict their influence on reaction rates.

Typical process equipment, key design issues, and perspectives on the current status of CFD for simulating reactive flows (Joshi and Ranade, 2003)

- 158 -

Appendix A

Typical equipment Design issues / applications

Monoliths, fixed beds, packed columns, trickle-bed reactors, filters, etc. Pressure drop, flow maldistribution and channelling, mixing and residence time distribution, heat/mass transfer coefficients, scale-up.

CURRENT STATUS Knowledge of physics

CFD codes / design applicability

Lumped models with or without isotropic porosity are used. Reasonable results for single phase flows; closure models suitable for multiphase flows through packed beds are not adequate. Commercial CFD codes/models allow simulations of single phase flows through porous media with isotropic or non-isotropic permeability and inertial resistance coefficients. Detailed modeling of a packed bed, which is very computational intensive, is being explored in recent studies. Capability of simulating multiphase flow through packed beds is almost nonexistent except via user defined functions based on empirically calibrated closure models. Used for designing equipment with single phase flow through a porous medium; use for designing multiphase flow through a packed bed is in the primitive stage.

PATH FORWARD Representing void space in packed beds: grid quality, computing resources.

Limitation of current models

Closures for multiphase flow through packed beds: interfacial area, drag. Simulation of flow regimes/transition and unsteady flows. Mixing and liquid dispersion: role of capillary forces/wetting.

Experiments

Development needed for overcoming the limitations

Physical models

Measurements of porosity distribution (with different length scales) within a packed bed (mean, axially averaged, standard deviation and so on) with quantification of the way of packing (Sharma et al., 2001). Pressure drop measurements along with local turbulence characteristics (Niven, 2002). Local phase distributions and velocity measurements for gas-liquid flow through well-characterized packed beds (Gladden et al., 2002). Detailed measurements of unsteady flow regimes (quantification of key spatial and temporal scales). Experimental data (global as well as local measurements) on liquid dispersion, bypass and channelling (Tsochatzidis et al., 2002).

Direct numerical simulations or VOF simulations are needed to guide the development of closures.

Models to simulate different packing possibilities (Spedding and Spencer, 1995). Interface closure models for estimating gassolid, gas-liquid and liquid-solid phases under different flow regimes (Jiang et al., 2002). Models to represent transient or low Re turbulent flow through complex geometry with severe curvature. Quantification of the role of gradients of porosity and capillary forces on liquid dispersion (Yin et al., 2000). Quantitative representation of the role of wetting and hysteresis on contact angle.

Table A-11.

Typical process equipment, key design issues, and perspectives on the current status of CFD for simulating flow through packed beds (Joshi and Ranade, 2003)

- 159 -

Appendix A

References Ambrosini, W., Forgione, N., Oriolo, F., (2002). Statistical characteristics of a water film falling down a flat plate at different inclinations and temperatures. International Journal of Multiphase Flow, 28, 1521 – 1540. Antunes, P., Gil, O., Bernardo-Gil, M.G., (2003). Supercritical fluid extraction of organochlorines from fish muscle with different sample preparation. Journal of Supercritical Fluids, 25, 135 - 142. Arunajatesan, V., Subramaniam, B., Hutchenson, K.W., Herkes, F.E., (2001). Fixed-bed hydrogenation of organic compounds in supercritical carbon dioxide. Chemical Engineering Science, 56, 1363 – 1369. Badalyan, A.G., Wilkinson, G.T., Chun, B.S., (1998). Extraction of Australian ginger root with carbon dioxide and ethanol entrainer. Journal of Supercritical Fluids, 13, 319 – 324. Baiker, A., (1999). Supercritical fluids in heterogeneous catalysis. Chemical Reviews, 99, 453 – 473. Barnabas, I.J., Dean, J.R., Tomlinson, W.R., Owen, S.P., (1995). Experimental design approach for the extraction of polycyclic aromatic hydrocarbons from soil using supercritical carbon dioxide. Analytical Chemistry, 67, 2064 – 2069. Baysal, T. and Starmans, D.A.J., (1999). Supercritical carbon dioxide extraction of carvone and limonene from caraway seed. Journal of Supercritical Fluids, 14, 225 – 234. Bertucco, A., Canu, P., Devetta, L., Zwahlen, A.G., (1997). Catalytic hydrogenation in supercritical CO2: kinetic measurements in a gradientless internal-recycle reactor. Industrial & Engineering Chemistry Research, 36, 2626 – 2633. Bombač, A., Žun, I., Filipič, B., Žumer, M., (1997). Gas-filled cavity structures and local void fraction distribution in aerated stirred vessel. AIChE Journal, 43, 2921 – 2931. Boselli, E. and Caboni, M.F., (2000). Supercritical carbon dioxide extraction of phospholipids from dried egg yolk without organic modifier. Journal of Supercritical Fluids, 19, 45 – 50. Brady, B.O., Kao, C.P.C., Dooley, K.M., Knopf, F.C., Gambrell, R.P., (1987). Supercritical extraction of toxic organics from soils. Industrial & Engineering Chemistry Research, 26, 261 – 268. Buchmann, M. and Mewes, D., (2000). Tomographic measurements of micro- and macromixing using the dual wavelength photometry. Chemical Engineering Journal, 77, 3 – 9. Buwa, V.V. and Ranade, V.V., (2002). Dynamics of gas-liquid flow in a rectangular bubble column: Experiments and single/multi-group CFD simulations. Chemical Engineering Science, 57, 4715 – 4736. Carlson, L.H.C., Machado, R.A.F., Spricigo, C.B., Pereira, L.K., Bolzan, A., (2001). Extraction of lemongrass essential oil with dense carbon dioxide. Journal of Supercritical Fluids, 21, 33 – 39.

- 160 -

Appendix A

Cassanello, M., Larachi, F., Kemoun, A., Al-Dahhan, M.H., Dudukovic, M.P., (2001). Inferring liquid chaotic dynamics in bubble columns. Chemical Engineering Science, 56, 6125 – 6134. Catchpole, O.J., Perry, N.B., da Silva, B.M.T., Grey, J.B., Smallfield, B.M., (2002). Supercritical extraction of herbs I: Saw Palmetto, St John's Wort, Kava Root, and Echinacea. Journal of Supercritical Fluids, 22, 129 – 138. Chiu, K.L., Cheng, Y.C., Chen, J.H., Chang, C.J., Yang, P.W., (2002). Supercritical fluids extraction of Ginkgo ginkgolides and flavonoids. Journal of Supercritical Fluids, 24, 77 – 87. Chouchi, D., Gourgouillon, D., Courel, M., Vital, J., Nunes da Ponte, M., (2001). The influence of phase behavior on reactions at supercritical conditions: The hydrogenation of α-pinene. Industrial & Engineering Chemistry Research, 40, 2551 – 2554. Churchill, S.W. and Zajic, S.C., (2002). Prediction of fully developed turbulent convection with minimal explicit empiricism. AIChE Journal, 48, 927 – 940. Ciofalo, M., Brucato, A., Grisafi, F., Torraca, N., (1996). Turbulent flow in closed and free-surface unbaffled tanks stirred by radial impellers. Chemical Engineering Science, 51, 3557 – 3573. Clark, M.C. and Subramaniam, B., (1996). 1-hexene isomerization on a Pt/γ-Al2O3 catalyst: The dramatic effects of feed peroxides on catalyst activity. Chemical Engineering Science, 51, 2369 – 2377. Clark, M.C. and Subramaniam, B., (1998). Extended alkylate production activity during fixedbed supercritical 1-butene/isobutane alkylation on solid acid catalysts using carbon dioxide as a diluent. Industrial & Engineering Chemistry Research, 37, 1243 – 1250. Collins, N.A., Debenedetti, P.G., Sundaresan, S., (1998). Disproportionation of toluene over ZSM-5 under near-critical conditions. AIChE Journal, 34, 1211 – 1214. Cooper-White, J.J., Fagan, J.E., Tirtaatmadja, V., Lester, D.R., Boger, D.V., (2002). Drop formation dynamics of constant low-viscosity, elastic fluids. Journal of Non-Newtonian Fluid Mechanics, 106, 29 – 59. Daood, H.G., Illés, V., Gnayfeed, M.H., Mészáros, B., Horváth, G., Biacs, A., (2002). Extraction of pungent spice paprika by supercritical carbon dioxide and subcritical propane. Journal of Supercritical Fluids, 23, 143 – 152. Dardas, Z., Süer, M.G., Ma, Y.H., Moser, W.R., (1996). High-temperature, high-pressure in situ reaction monitoring of heterogeneous catalytic processes under supercritical conditions by CIR-FTIR. Journal of Catalysis, 159, 204 – 211. Daukšas, E., Venskutonis, P.R., Sivik, B., (1999). Supercritical CO2 extraction of the main constituents of lovage (Levisticum officinale Koch.) essential oil in model systems and overground botanical parts of the plant. Journal of Supercritical Fluids, 15, 51 – 62. Daukšas, E., Venskutonis, P.R., Sivik, B., (2002). Supercritical fluid extraction of borage (Borago officinalis L.) seeds with pure CO2 and its mixture with caprylic acid methyl ester. Journal of Supercritical Fluids, 22, 211 – 219. de Azevedo, A.B.A., Kopcak, U., Mohamed, R.S., (2003). Extraction of fat from fermented Cupuaçu seeds with supercritical solvents. Journal of Supercritical Fluids, 27, 223 – 237.

- 161 -

Appendix A

del Valle, J.M., Jiménez, J.C., de la Fuente, J.C., (2003). Extraction kinetics of pre-pelletized Jalapeño peppers with supercritical CO2. Journal of Supercritical Fluids, 25, 33 – 44. Desai, R.B. and Joshi, J.B., (1996). Flow pattern in the vicinity of gas-liquid interface in a stirred cell. Journal of Chemical Engineering of Japan, 29, 568 – 576. Doneanu, C. and Anitescu, G., (1998). Supercritical carbon dioxide extraction of Angelica archangelica L. root oil. Journal of Supercritical Fluids, 12, 59 – 67. Dooley, K.M. and Knopf, F.C., (1987). Oxidation catalysis in a supercritical fluid medium. Industrial and Engineering Chemistry Research, 26, 1910 – 1916. Dunford, N.T., Goto, M., Temelli, J., (1998). Modeling of oil extraction with supercritical CO2 from Atlantic mackerel (Scomber scombrus) at different moisture contents. Journal of Supercritical Fluids, 13, 303 – 309. Esquível, M.M., Ribeiro, M.A., Bernardo-Gil, M.G., (1999). Supercritical extraction of savory oil: study of antioxidant activity and extract characterization. Journal of Supercritical Fluids, 14, 129 – 138. Eychenne, V., Sáiz, S., Trabelsi, F., Recasens, F., (2001). Near-critical solvent extraction of wool with modified carbon dioxide — experimental results. Journal of Supercritical Fluids, 21, 23 – 31. Fan, L., Nakamura, I., Ishida, S., Fujimoto, K., (1997a). Supercritical-phase alkylation reaction on solid acid catalysts: mechanistic study and catalyst development. Industrial & Engineering Chemistry Research, 36, 1458 – 1463. Fan, L., Watanabe, T., Fujimoto, K., (1997b). Reaction phase effect on tertiary butyl alcohol synthesis by air oxidation of isobutane. Applied Catalysis A: General, 158, L41 - L46. Ferreira, S.R.S., Nikolov, Z.L., Doraiswamy, L.K., Meireles, M.A.A., Petenate, A.J., (1999). Supercritical fluid extraction of black pepper (Piper nigrun L.) essential oil. Journal of Supercritical Fluids, 14, 235 – 245. Fourcade, E., Wadley, R., Hoefsloot, H.C.J., Green, A., Iedema, P.D., (2001). CFD calculation of laminar striation thinning in static mixer reactors. Chemical Engineering Science, 56, 6729 – 6741. Fukai, J., Shiiba, Y., Yamamoto, T., Miyatake, O., Pouliakos, D., Megaridis, C., Zhao, Z., (1995). Wetting effects on the spreading of a liquid droplet colliding with a flat surface: Experiment and modeling. Physics of Fluids, 7, 236 – 247. Gladden, L.F., Mantle, M.D., Sederman, A.J., Yuen, E.H.L., (2002). Magnetic resonance imaging of single- and two-phase flow in fixed-bed reactors. Applied Magnetic Resonance, 22, 201 – 212. Gouesbet, G. and Berlemont, A., (1999). Eulerian and Lagrangian approaches for predicting the behaviour of discrete particles in turbulent flows. Progress in Energy and Combustion Science, 25, 133 – 159. Gulawani, S.S., (2002). Hydrodynamics at interfaces and theories of mass transfer. Ph.D. Thesis. University of Mumbai, Mumbai.

- 162 -

Appendix A

Harris, A.T., Davidson, J.F., Thorpe, R.B., (2002). The prediction of particle cluster properties in the near wall region of a vertical riser (200157). Powder Technology, 127, 128 – 143. Hitzler, M.G., Poliakoff, M., (1997). Continuous hydrogenation of organic compounds in supercritical fluids. Chemical Communications, 1667 – 1668. Hitzler, M.G., Smail, F.R., Poliakoff, M., Ross, S.K., (1998). Friedel–Crafts alkylation in supercritical fluids: continuous, selective and clean. Chemical Communications, 359 – 360. Illés, V., Szalai, O., Then, M., Daood, H., Perneczki, S., (1997). Extraction of hiprose fruit by supercritical CO2 and propane. Journal of Supercritical Fluids, 10, 209 – 218. Illés, V., Daood, H.G., Perneczki, S., Szokonya, L., Then, J., (2000). Extraction of coriander seed oil by CO2 and propane at super- and subcritical conditions. Journal of Supercritical Fluids, 17, 177 – 186. Izquierdo, A., Tena, M.T., de Castro, M.L.D., Valcarcel, M., (1996). Supercritical fluid extraction of carbamate pesticides from soils and cereals. Cromatographia, 42, 206 – 212. Jana, S.C., Tjahjadi, M., Ottino, J.M., (1994). Chaotic mixing of viscous fluids by periodic changes in geometry: Baffled cavity flow. AIChE Journal, 40, 1769 – 1781. Jessop, P.G. and Leitner, W., (1999). Chemical synthesis using supercritical fluids. Wiley - VCH, Weinheim. Jiang, Y., Khadilkar, M.R., Al-Dahhan, M.H., Dudukovic, M.P., (2002). CFD of multiphase flow in packed-bed reactors: I. k-Fluid modeling issues. AIChE Journal, 48, 701 – 715. Joshi, J.B., Sawant, S.B., Patwardhan, A.W., Patil, D.J., Kshatriya, S.S., Nere, N.K., (2001). Relation between flow pattern and de-activation of enzymes in stirred reactors. Chemical Engineering Science, 56, 443 – 452. Joshi, J. and Ranade, V., (2003). Computational fluid dynamics for designing process equipment: expectations, current status, and path forward. Industrial & Engineering Chemistry Research, 42, 1115 - 1128. Karcz, J., (1999). Studies of local heat transfer in a gas-liquid system agitated by double disc turbines in a slender vessel. Chemical Engineering Journal, 72, 217 – 227. King, J.W., Holliday, R.L., List, G.R., Snyder, J.M., (2001). Hydrogenation of vegetable oils using mixtures of supercritical carbon dioxide and hydrogen. JAOCS, Journal of the American Oil Chemists' Society, 78, 107 – 113. Kiriamiti, H.K., Camy, S., Gourdon, C., Condoret, J.S., (2003). Pyrethrin extraction from pyrethrum flowers using carbon dioxide. Journal of Supercritical Fluids, 26, 193 – 200. Kresta, S.M. and Wood, P.E., (1993). The flow field produced by a pitched blade turbine: Characterization of the turbulence and estimation of the dissipation rate. Chemical Engineering Science, 48, 1761 – 1774. Liou, T.-M., Chen, S.-H., Shih, K.-C., (2002). Numerical simulation of turbulent flow field and heat transfer in a two-dimensional channel with periodic slit ribs. International Journal of Heat and Mass Transfer, 45, 4493 – 4505.

- 163 -

Appendix A

Lojková, L., Sedláková, J., Kubáň, V., (2005). A two-step supercritical fluid extraction of polycyclic aromatic hydrocarbons from roadside soil samples. Journal of Separation Science, 28, 2067 – 2075. López-Mesas, M., Christoe, J., Carrillo, F., Crespi, M., (2005). Supercritical fluid extraction with cosolvents of wool wax from wool scour wastes. Journal of Supercritical Fluids, 35, 235 – 239. Lu, W.-M. and Ju, S.-J., (1989). Cavity configuration, flooding and pumping capacity of disctype turbines in aerated stirred tanks. Chemical Engineering Science, 44, 333 – 342. Macher, M.-B., Holmqvist, A., (2001). Hydrogenation of palm oil in near-critical and supercritical propane. European Journal of Lipid Science and Technology, 103, 81 – 84. Marchisio, D.L., Fox, R.O., Barresi, A.A., Garbero, M., Baldi, G., (2001). On the simulation of turbulent precipitation in a tubular reactor via computational fluid dynamics (CFD). Chemical Engineering Research and Design, 79, 998 – 1004. Marsal, A., Celma, P.J., Cot, J., Cequier, M., (2000a). Supercritical CO2 extraction as a clean degreasing process in the leather industry. Journal of Supercritical Fluids, 16, 217 – 223. Marsal, A., Celma, P.J., Cot, J., Cequier, M., (2000b). Application of the supercritical CO2 extraction technology on the recovery of natural fat from the sheepskin degreasing process. Journal of Supercritical Fluids, 18, 65 – 72. McCoy, B.J. and Subramaniam, B., (1995). Continuous-mixture kinetics of coke formation from olefinic oligomers. AIChE Journal, 41, 317 – 323. McDaniel, L.H., Ashraf-Khorassani, M., Taylor, L.T., (2001). Supercritical fluid extraction of wood pulp with analysis by capillary gas chromatography-mass spectrometry. Journal of Supercritical Fluids, 19, 275 – 286. Mira, B., Blasco, M., Berna, A., Subirats, S., (1999). Supercritical CO2 extraction of essential oil from orange peel. Effect of operation conditions on the extract composition. Journal of Supercritical Fluids, 14, 95 – 104. Montanari, L., Fantozzi, P., Snyder, J.M., King, J.W., (1999). Selective extraction of phospholipids from soybeans with supercritical carbon dioxide and ethanol. Journal of Supercritical Fluids, 14, 87 – 93. Monteiro, A.R., Meireles, M.A.A., Marques, M.O.M, Petenate, A.J., (1997). Extraction of the soluble material from the shells of the bacuri fruit (Platonia insignis Mart) with pressurized CO2 and other solvents. Journal of Supercritical Fluids, 11, 91 – 102. Niu, F. and Hofmann, H., (1995). Deactivation kinetics and modelling of coke removal under supercritical conditions for the example of ethylbenzene disproportionation. Chemical Engineering & Technology, 18, 278 – 283. Niu, F. and Hofmann, H., (1997). Investigation of coke extraction from zeolite-HY under supercritical and near-critical conditions. Canadian Journal of Chemical Engineering, 75, 346 – 352. Niven, R.K., (2002). Physical insight into the Ergun and Wen and Yu equations for fluid flow in packed and fluidised beds. Chemical Engineering Science, 57, 527 – 534.

- 164 -

Appendix A

Palikaras, A., Yakinthos, K., Goulas, A., (2002). Transition on a flat plate with a semi-circular leading edge under uniform and positive shear free-stream flow. International Journal of Heat and Fluid Flow, 23, 455 – 470. Pasquini, D., Pimenta, M.T.B., Ferreira, L.H., Curvelo, A.A. d. S., (2005). Extraction of lignin from sugar cane bagasse and Pinus taeda wood chips using ethanol–water mixtures and carbon dioxide at high pressures. Journal of Supercritical Fluids, 36, 31 – 39. Povh, N.P., Marques, M.O.M., Meireles, M.A.A., (2001). Supercritical CO2 extraction of essential oil and oleoresin from chamomile (Chamomilla recutita [L.] Rauschert). Journal of Supercritical Fluids, 21, 245 – 256. Ramírez, E., Zgarni, S., Larrayoz, M.A., Recasens, F., (2002). Short compilation of published reaction rate data for catalytic hydrogenations in supercritical fluids. Chemical Engineering and Technology: Engineering in Life Science, 25, 257 – 264. Ramírez, E., Recasens, F., Fernández, M., Larrayoz, M.A., (2004). Sunflower oil hydrogenation on Pd/C in SC propane in a continuous recycle reactor. AIChE Journal, 50, 1545 – 1555. Ramírez, E., (2005). Contribution to the study of heterogeneous catalytic reactions in SCFs: hydrogenation of sunflower oil in Pd catalysts at single-phase conditions. Ph.D. Thesis. Universitat Politècnica de Catalunya, Barcelona. Ribeiro, M.A., Bernardo-Gil, M.G., Esquivel, M.M., (2001). Melissa officinalis, L.: study of antioxidant activity in supercritical residues. Journal of Supercritical Fluids, 21, 51 – 60. Rónyai, E., Simándi, B., Tömösközi, S., Deák, A., Vigh, L., Weinbrenner, Z., (1998). Supercritical fluid extraction of corn germ with carbon dioxide–ethyl alcohol mixture. Journal of Supercritical Fluids, 14, 75 – 81. Rothstein, J.P. and McKinley, G.H., (2001). The axisymmetric contraction-expansion: The role of extensional rheology on vortex growth dynamics and the enhanced pressure drop. Journal of Non-Newtonian Fluid Mechanics, 98, 33 – 63. Saldaña, M.D.A., Zetzl, C., Mohamed, R.S., Brunner, G., (2002). Decaffeination of guaraná seeds in a microextraction column using water-saturated CO2. Journal of Supercritical Fluids, 22, 119 – 127. Sankaranarayanan, K. and Sundaresan, S., (2002). Lift force in bubbly suspensions. Chemical Engineering Science, 57, 3521 – 3542. Shamsipur, M., Ghiasvand, A.R., Yamini, Y., (2001). Extraction of uranium from solid matrices using modified supercritical fluid CO2. Journal of Supercritical Fluids, 20, 163 – 169. Sharma, S., Mantle, M.D., Gladden, L.F., Winterbottom, J.M., (2001). Determination of bed voidage using water substitution and 3D magnetic resonance imaging, bed density and pressure drop in packed-bed reactors. Chemical Engineering Science, 56, 587 – 595. Simándi, B., Kristo, Sz.T., Kéry, Á., Selmeczi, L.K., Kmecz, I., Kemény, S., (2002). Supercritical fluid extraction of dandelion leaves. Journal of Supercritical Fluids, 23, 135 – 142. Snavely, K. and Subramaniam, B., (1997). On-line gas chromatographic analysis of FischerTropsch synthesis products formed in a supercritical reaction medium. Industrial & Engineering Chemistry Research, 36, 4413 – 4420.

- 165 -

Appendix A

Son, S.M. and Singh, R.K., (2002). Turbulence modeling and verification for aseptically processed soybean milk under turbulent flow conditions. Journal of Food Engineering, 52, 177 – 184. Spedding, P.L. and Spencer, R.M., (1995). Simulation of packing density and liquid flow in fixed beds. Computers and Chemical Engineering, 19, 43 – 73. Speziale, C.G., (1998). Turbulence modeling for time-dependent RANS and VLES: A review. AIAA Journal, 36, 173 – 184. Spricigo, C.B., Pinto, L.T., Bolzan, A., Novais, A.F., (1999). Extraction of essential oil and lipids from nutmeg by liquid carbon dioxide. Journal of Supercritical Fluids, 15, 253 – 259. Subra, P., Castellani, S., Jestin, P., Aoufi, A., (1998). Extraction of β − carotene with supercritical fluids: Experiments and modeling. Journal of Supercritical Fluids, 12, 261 – 269. Suga, K., (2003). Predicting turbulence and heat transfer in 3-D curved ducts by near-wall second moment closures. International Journal of Heat and Mass Transfer, 46, 161 – 173. Sun, X., Smith, T.R., Kim, S., Ishii, M., Uhle, J., (2003). Interfacial area of bubbly flow in a relatively large diameter pipe. Experimental Thermal and Fluid Science, 27, 97 – 109. Tacke, T., Wieland, S., Panster, P., (2003). Selective and complete hydrogenation of vegetable oils and free fatty acids in supercritical fluids. In “Green chemistry using liquid and supercritical carbon dioxide”, ed. J. M. De Simone. Oxford University Press, Oxford. Ted Mao, D.C., Kuhn, S., Tran, H. Spread and rebound of liquid droplets upon impact on flat surfaces. AIChE Journal, 43, 2169 – 2179. Tiltscher, H., Wolf, H., Schelschshorn, J., (1984). Utilization of supercritical fluid solvent-effects in heterogeneous catalysis. Berichte der Bunsengesellschaft/Physical Chemistry Chemical Physics, 88, 897 – 900. Tiltscher, H. and Hofmann, H., (1986). Trends in high pressure chemical reaction engineering. Chemical Engineering Science, 42, 959 – 977. Tomiyama, A., Tamai, H., Zun, I., Hosokawa, S., (2002). Transverse migration of single bubbles in simple shear flows. Chemical Engineering Science, 57, 1849 – 1858. Tsochatzidis, N.A., Karabelas, A.J., Giakoumakis, D., Huff, G.A., (2002). An investigation of liquid maldistribution in trickle beds. Chemical Engineering Science, 57, 3543 – 3555. van der Bosch, R.H., (1998). The role of clusters in gas-solid reactors: An experimental study. Ph.D. Thesis. Twente University, Enschede. van den Hark, S., Härröd, M., Møller, P., (1999). Hydrogenation of fatty acid methyl esters to fatty alcohols at supercritical conditions. JAOCS, Journal of the American Oil Chemists' Society, 76, 1363 – 1370. van den Hark, S., Härröd, M., Møller, P., (2000). Production of fatty alcohols by heterogeneous catalysis at supercritical single-phase conditions. Studies in Surface Science and Catalysis, 130 A, 539 – 544.

- 166 -

Appendix A

van den Hark, S., Härröd, M., (2001). Fixed-bed hydrogenation at supercritical conditions to form fatty alcohols: The dramatic effects caused by phase transitions in the reactor. Industrial & Engineering Chemistry Research, 40, 5052 – 5057. van Ruymbeke, E., Keunings, R., Stéphenne, V., Hagenaars, A., Bailly, C., (2002). Evaluation of reptation models for predicting the linear viscoelastic properties of entangled linear polymers. Macromolecules, 35, 2689 – 2699. Vaquero, E., Beltrán, S., Sanz, M.T., (2006). Extraction of fat from pigskin with supercritical carbon dioxide. Journal of Supercritical Fluids, in Press. Vedaraman, N., Srinivasakannan, C., Brunner, G., Ramabrahmam, B.V., Rao, P.G., (2005). Experimental and modeling studies on extraction of cholesterol from cow brain using supercritical carbon dioxide. Journal of Supercritical Fluids, 34, 27 – 34. Verschuren, I.L.M., Wijers, J.G., Keurentjes, J.T.F., (2002). Mean concentrations and concentration fluctuations in a stirred-tank reactor. AIChE Journal, 48, 1390 – 1400. Vieville, C., Mouloungui, Z., Gaset, A., (1993). Esterification of oleic acid by methanol catalyzed by p-toluenesulfonic acid and the cation-exchange resins K2411 and K1481 in supercritical carbon dioxide. Industrial & Engineering Chemistry Research, 32, 2065 – 2068. Vogel, J.C., (1984). Heat transfer and fluid mechanics measurements in the turbulent reattaching flow behind a backward-facing flow step. Ph.D. Thesis. Stanford University, Stanford. Yang, Y., Hawthorne, S.B., Miller, D.J., (1995). Comparison of sorbent and solvent trapping after supercritical fluid extraction of volatile petroleum hydrocarbons from soil. Journal of Chromatography A, 699, 265 – 276. Yin, F.H., Sun, C.G., Afacan, A., Nandakumar, K., Chuang, K.T., (2000). CFD modeling of masstransfer processes in randomly packed distillation columns. Industrial & Engineering Chemistry Research, 39, 1369 – 1380. Yokota, K., Hanakata, Y., Fujimoto, K., (1990). Supercritical phase Fischer-Tropsch synthesis. Chemical Engineering Science, 45, 2743 – 2750. Yokota, K. and Fujimoto, K., (1991a). Supercritical-phase Fischer-Tropsch synthesis reaction. 2. The effective diffusion of reactant and products in the supercritical-phase reaction. Industrial & Engineering Chemistry Research, 30, 95 – 100. Yokota, K., Hanakata, Y., Fujimoto, K., (1991b). Supercritical phase Fischer-Tropsch synthesis reaction. 3. Extraction capability of supercritical fluids. Fuel, 70, 989 – 994. Zancan, K.C., Marques, M.O.M., Petenate, A.J., Meireles, M.A.A., (2002). Extraction of ginger (Zingiber officinale Roscoe) oleoresin with CO2 and co-solvents: a study of the antioxidant action of the extracts. Journal of Supercritical Fluids, 24, 57 – 76.

- 167 -

APPENDIX B

TURBULENCE MODELS

Transport equations, turbulent quantities and convective heat and mass transfer modeling

Michael Graf, « Turbulence »

Appendix B

B-1. The Spalart-Allmaras model B-1.1. Transport equation In turbulence models that employ the Boussinesq approach, the central issue is how the eddy viscosity is computed. The model proposed by Spalart and Allmaras (1992) solves a transport equation for a quantity that is a modified form of the turbulent kinematic viscosity. The transported variable in the Spalart-Allmaras model, ν~ , is identical to the turbulent kinematic viscosity except in the near-wall (viscous-affected) region. The transport equation for ν~ is 2 ⎛ ∂ν~ ⎞ ⎤ ∂ ∂ 1 ⎡⎢ ∂ ⎧⎪ ∂ν~ ⎫⎪ ~ ~ ~ ⎟ ⎥ − Yν + Sν~ (ρ ν ) + (ρ ν u i ) = Gν + ⎨(µ + ρ ν ) ⎬ + C b 2 ρ ⎜⎜ ⎟ ∂t ∂x i ∂x j ⎪⎭ ∂ x σ ν~ ⎢ ∂x j ⎪⎩ ⎝ j ⎠ ⎥⎦ ⎣

[B-1] where Gν is the production of turbulent viscosity and Yν is the destruction of turbulent viscosity that occurs in the near-wall region due to wall blocking and viscous damping. σ ν~ and C b 2 are constants and ν is the molecular kinematic viscosity. Sν~ is a user-defined source term. Note that since the turbulence kinetic energy is not calculated in the Spalart-Allmaras model, the last term in Eq. [2.1-8] is ignored when estimating the Reynolds stresses.

B-1.2. Modeling the turbulent viscosity The turbulent viscosity, µ t , is computed from

µ t = ρν~f v1

[B-2]

where the viscous damping function, f v1 , is given by

f v1 =

Χ≡

Χ3 Χ 3 + C v31

[B-3]

ν~ ν

[B-4]

B-1.3. Modeling the turbulent production The production term, Gν , is modeled as

~ Gν = C b1 ρ S ν~

[B-5]

ν~ ~ Θ ≡ Θ + 2 2 f v2 κ d

[B-6]

- 171 -

Appendix B

f v2 = 1 −

Χ 1 + Χ f v1

[B-7]

C b1 and κ are constants, d is the distance from the wall, and Θ is a scalar measure of the

deformation tensor. By default, Θ is based on the magnitude of the vorticity:

Θ ≡ 2Ω ij Ω ij

[B-8]

where Ω ij is the mean rate-of-rotation tensor and is defined by

Ω ij =

1 ⎛⎜ ∂u i ∂u j − 2 ⎜⎝ ∂x j ∂xi

⎞ ⎟ ⎟ ⎠

[B-9]

B-1.4. Modeling the turbulent destruction The destruction term is modeled as

⎛ ν~ ⎞ Yν = C w1 ρ f w ⎜ ⎟ ⎝d ⎠

2

[B-10]

⎡ 1+ C6 ⎤ f w = g ⎢ 6 w36 ⎥ ⎣ g + C w3 ⎦

(

g = r + C w2 r 6 − r

1

6

[B-11]

)

[B-12]

ν~ r≡ ~ 2 2 Θκ d

[B-13]

~

C w1 , C w2 and C w3 are constants, and Θ is given by Eq. [B-6].

B-1.5. Model constants The model constants have the following default values (Spalart and Allmaras, 1992):

C b1 = 0.1355; C w1 =

C b1

κ2

+

C b 2 = 0.622;

1 + Cv 2

σ ν~

;

2 3

σ ν~ = ;

C w 2 = 0.3;

C v1 = 7.1;

C w3 = 2.0;

- 172 -

κ = 0.4187;

[B-14]

Appendix B

B-1.6. Wall boundary conditions At walls, the modified turbulent kinematic viscosity, ν~ , is set to zero. When the mesh is fine enough to resolve the laminar sublayer, the wall shear stress is obtained from the laminar stress-strain relationship:

u ρ ut d = ut µ

[B-15]

If the mesh is too coarse to resolve the laminar sublayer, it is assumed that the centroid of the wall-adjacent cell falls within the logarithmic region of the boundary layer, and the law-of-thewall is employed:

⎛ ρ ut d ⎞ 1 u ⎟⎟ = ln 9.793⎜⎜ u t 0.4187 ⎝ µ ⎠

[B-16]

where u is the velocity parallel to the wall and

u t is the shear velocity.

B-1.7. Convective heat and mass transfer modeling Turbulent heat transport is modeled using the concept of Reynolds' analogy to turbulent momentum transfer. The "modeled'' energy equation is thus given by the following:

⎡⎛ C µ ∂ (ρE ) + ∂ [ui (ρE + p )] = ∂ ⎢⎜⎜ k + p t ∂t ∂xi ∂x j ⎢⎣⎝ Prt where k is the thermal conductivity,

⎤ ⎞ ∂T ⎟⎟ + u i (τ ij )eff ⎥ + S h ⎥⎦ ⎠ ∂x j

[B-17]

E is the total energy, and (τ ij )eff is the deviatoric stress

tensor, defined as

(τ )

ij eff

⎛ ∂u j ∂u i = µ eff ⎜ + ⎜ ∂x ⎝ i ∂x j

( )eff

The term involving τ ij

⎞ 2 ⎟ − µ eff ∂u i δ ij ⎟ 3 ∂xi ⎠

[B-18]

represents the viscous heating. The default value of the turbulent

Prandtl number is 0.85. Turbulent mass transfer is treated similarly, with a default turbulent Schmidt number of 0.7. Wall boundary conditions for scalar transport are handled analogously to momentum, using the appropriate "law-of-the-wall''.

B-2. The κ − ε family models This section presents the standard, RNG, and realizable κ − ε models. All three models have similar forms, with transport equations for κ and ε . The major differences in the models are as follows:

- 173 -

Appendix B

•

the method of calculating turbulent viscosity

•

the turbulent Prandtl numbers governing the turbulent diffusion of κ and ε .

•

the generation and destruction terms in the equation.

The transport equations, methods of calculating turbulent viscosity, and model constants are presented separately for each model. The features that are essentially common to all models follow, including turbulent production, generation due to buoyancy, accounting for the effects of compressibility, and modeling heat and mass transfer.

B-2.1. The standard κ − ε model The standard κ − ε model (Launder and Spalding, 1972) is a semi-empirical model based on model transport equations for the turbulence kinetic energy ( κ ) and its dissipation rate ( ε ). The model transport equation for κ is derived from the exact equation, while the model transport equation for ε was obtained using physical reasoning and bears little resemblance to its mathematically exact counterpart. In the derivation of the κ − ε model, it was assumed that the flow is fully turbulent, and the effects of molecular viscosity are negligible. The standard κ − ε model is therefore valid only for fully turbulent flows.

B-2.1.1. Transport equations The turbulence kinetic energy, κ , and its rate of dissipation, ε , are obtained from the following transport equations:

⎡⎛ ∂ (ρ κ ) + ∂ (ρ κ u i ) = ∂ ⎢⎜⎜ µ + µ t (Prt )κ ∂t ∂xi ∂x j ⎣⎢⎝ ⎡⎛ ∂ (ρε ) + ∂ (ρε u i ) = ∂ ⎢⎜⎜ µ + µ t (Prt )ε ∂t ∂xi ∂x j ⎢⎣⎝

⎞ ∂κ ⎤ ⎟ ⎟ ∂x ⎥ + Gκ − ρε + S κ ⎠ j ⎦⎥

[B-19]

⎞ ∂ε ⎤ ε ε2 ⎟ ( ) + C G + C G − C + Sε ρ 1ε 3ε b 2ε k ⎟ ∂x ⎥ κ κ ⎥ j ⎠ ⎦ [B-20]

In these equations, Gκ represents the generation of turbulence kinetic energy due to the mean velocity gradients, Gb is the generation of turbulence kinetic energy due to buoyancy, Ym represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate, C1ε , C 2ε , and C 3ε are constants, and S κ and S ε are user-defined source terms.

B-2.1.2. Modeling the turbulent viscosity The turbulent (or eddy) viscosity , µ t , is computed by combining

- 174 -

κ

and ε as follows:

Appendix B

µt = ρ Cµ

κ2 ε

[B-21]

where C µ is a constant.

B-2.1.3. Model constants The model constants have the following default values (Launder and Spalding, 1972):

C1ε = 1.44;

C 2ε = 1.92;

C µ = 0.09;

(Prt )κ

= 1.0;

(Prt )ε

= 1.3;

[B-22]

These default values have been determined from experiments with air and water for fundamental turbulent shear flows including homogeneous shear flows and decaying isotropic grid turbulence. They have been found to work fairly well for a wide range of wall-bounded and free shear flows.

B-2.2. The RNG κ − ε model The RNG-based κ − ε turbulence model (Choudhury et al., 1993) is derived from the instantaneous Navier-Stokes equations, using a mathematical technique called “renormalization group” (RNG) methods. The analytical derivation results in a model with constants different from those in the standard κ − ε model, and additional terms and functions in the transport equations for κ and ε .

B-2.2.1. Transport equations The RNG κ − ε model has a similar form to the standard κ − ε model:

⎛ ⎞ ∂ (ρκ ) + ∂ (ρκ ui ) = ∂ ⎜⎜α κ µ eff ∂κ ⎟⎟ + Gκ + Gb − ρε − YM + S κ ∂t ∂xi ∂x j ⎝ ∂x j ⎠

[B-23]

2 ⎛ ⎞ ∂ (ρε ) + ∂ (ρε u i ) = ∂ ⎜⎜ α ε µ eff ∂ε ⎟⎟ + C1ε ε (Gκ + C3ε Gb ) − C 2ε ρ ε − Rε + S ε ∂t ∂xi ∂x j ⎝ ∂x j ⎠ κ κ

[B-24] In these equations, Gκ represents the generation of turbulence kinetic energy due to the mean velocity gradients, Gb is the generation of turbulence kinetic energy due to buoyancy, and Ym represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. The quantities α κ and α ε are the inverse effective Prandtl numbers for κ and ε , respectively. S κ and S ε are user-defined source terms.

- 175 -

Appendix B

B-2.2.2. Modeling the effective viscosity The scale elimination procedure in RNG theory results in a differential equation for turbulent viscosity:

⎛ ρ 2κ d⎜ ⎜ εµ ⎝

νˆ =

⎞ νˆ ⎟ = 1.72 dνˆ ⎟ νˆ 3 − 1 + C v ⎠

[B-25]

µ eff µ

[B-26]

Cν ≈ 100

[B-27]

Eq. [B-25] is integrated to obtain an accurate description of how the effective turbulent transport varies with the effective Reynolds number (or eddy scale), allowing the model to better handle low-Reynolds-number and near-wall flows . In the high-Reynolds-number limit, Equation [B-25] gives

µt = ρ Cµ

κ2 ε

[B-28]

with C µ = 0.0845, derived using RNG theory. It is interesting to note that this value of C µ is very close to the empirically-determined value of 0.09 used in the standard κ − ε model.

B-2.2.3. Calculating the inverse effective Prandtl numbers The inverse effective Prandtl numbers, α κ and α ε , are computed using the following formula derived analytically by the RNG theory:

α − 1.3929 α 0 − 1.3929

0.6321

α + 2.3929 α 0 + 2.3929

0.3679

=

µ mol µ eff

[B-29]

where α 0 = 1.0. In the high-Reynolds-number limit ( µ mol / µ eff η 0 ), however, the Rε term makes a negative contribution, making the value of C 2*ε less than C 2ε . In comparison with the standard κ − ε model, the smaller destruction of ε augments ε , reducing κ and, eventually, the effective viscosity. As a result, in rapidly strained flows, the RNG model yields a lower turbulent viscosity than the standard κ − ε model. Thus, the RNG model is more responsive to the effects of rapid strain and streamline curvature than the standard κ − ε model, which explains the superior performance of the RNG model for certain classes of flows.

B-2.2.5. Model constants The model constants have values derived analytically by the RNG theory. These values are:

C1ε = 1.42;

C 2ε = 1.68;

[B-34]

B-2.3. The realizable κ − ε model In addition to the standard and RNG-based κ − ε models, Shih et al. (1995) proposed the socalled realizable κ − ε model. The term “realizable” means that the model satisfies certain mathematical constraints on the normal stresses, consistent with the physics of turbulent flows. To understand this, consider combining the Boussinesq relationship (Eq. [2.1-8]) and the eddy viscosity definition (Eq. [B-21]) to obtain the following expression for the normal Reynolds stress in an incompressible strained mean flow:

∂U 2 u 2 = κ − 2ν t ∂x 3

[B-35]

- 177 -

Appendix B

Using Eq. [B-21] for ν t ≡ µ t / ρ , one obtains the result that the normal stress, u 2 , which by definition is a positive quantity, becomes negative, i.e., “non-realizable”, when the strain is large enough to satisfy

κ ∂U 1 ≥ ≈ 3 .7 ε ∂x 3C µ

[B-36]

Similarly, it can also be shown that the Schwarz inequality for shear stresses ( u α u β

2

≤ u α2 u β2 ; no

summation over α and β ) can be violated when the mean strain rate is large. The most straightforward way to ensure the realizability (positivity of normal stresses and Schwarz inequality for shear stresses) is to make C µ variable by sensitizing it to the mean flow (mean deformation) and the turbulence ( κ , ε ). The notion of variable C µ is suggested by many modelers including Reynolds (1987), and is well substantiated by experimental evidence. For example, C µ is found to be around 0.09 in the inertial sublayer of equilibrium boundary layers, and 0.05 in a strong homogeneous shear flow. Another weakness of the standard κ − ε model or other traditional κ − ε models lies with the modeled equation for the dissipation rate. The well-known round-jet anomaly (named based on the finding that the spreading rate in planar jets is predicted reasonably well, but prediction of the spreading rate for axisymmetric jets is unexpectedly poor) is considered to be mainly due to the modeled dissipation equation. The realizable κ − ε model was intended to address these deficiencies of traditional κ − ε models by adopting the following: •

a new eddy-viscosity formula involving a variable originally proposed by Reynolds (1987).

•

a new model equation for dissipation based on the dynamic equation of the meansquare vorticity fluctuation.

B-2.3.1. Transport equations The modeled transport equations for κ and ε in the realizable κ − ε model are

⎡⎛ ∂ (ρκ ) + ∂ (ρκ u j ) = ∂ ⎢⎜⎜ µ + µ t (Prt )k ∂t ∂x j ∂x j ⎣⎢⎝

⎞ ∂κ ⎤ ⎟ ⎟ ∂x ⎥ + Gκ + Gb − ρε − YM + S κ ⎠ j ⎦⎥

⎡⎛ ∂ (ρε ) + ∂ (ρε u j ) = ∂ ⎢⎜⎜ µ + µ t (Prt )ε ∂t ∂x j ∂x j ⎢⎣⎝

⎞ ∂ε ⎤ ε2 ε ⎟ ρ ε ρ + C V − C + C C 3ε Gb + S ε ⎥ 1 2 1 ε ⎟ ∂x κ κ νε + ⎥ j ⎠ ⎦

[B-37]

[B-38]

⎡ η ⎤ C1 = max ⎢0.43, ; η + 5 ⎥⎦ ⎣

η =V

κ ; ε

V = 2VijVij

- 178 -

[B-39]

Appendix B

In these equations, Gκ represents the generation of turbulence kinetic energy due to the mean velocity gradients, Gb is the generation of turbulence kinetic energy due to buoyancy, and Ym represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. S κ and S ε are user-defined source terms. Note that the κ equation (Eq. [B-37]) is the same as that in the standard κ − ε model (Eq. [B19]) and the RNG κ − ε model (Eq. [B-23]), except for the model constants. However, the form of the ε equation is quite different from those in the standard and RNG-based κ − ε models (Eqs. [B-20] and [B-24]). One of the noteworthy features is that the production term in the equation (the second term on the right-hand side of Eq. [B-38]) does not involve the production of κ ; i.e., it does not contain the same Gκ term as the other κ − ε models. It is believed that the present form better represents the spectral energy transfer. Another desirable feature is that the destruction term (the next to last term on the right-hand side of Eq. [B-38]) does not have any singularity; i.e., its denominator never vanishes, even if κ vanishes or becomes smaller than zero. This feature is contrasted with traditional κ − ε models, which have a singularity due to κ in the denominator.

B-2.3.2. Modeling the turbulent viscosity As in other κ − ε models, the turbulent viscosity is computed from:

κ2 µ t = ρC µ ε

[B-40]

The difference between the realizable κ − ε model and the standard and RNG κ − ε models is that C µ is no longer constant. It is computed from

Cµ =

1

[B-41]

κU * Ao + As ε

~ ~ U * ≡ VijVij + Ω ij Ω ij

[B-42]

~ Ω ij = Ω ij − 2ε ijk ω k

[B-43]

Ω ij = Ω ij − ε ijk ω k

[B-44]

where Ω ij is the mean rate-of-rotation tensor viewed in a rotating reference frame with the angular velocity ω k . The model constants A0 and As are given by

As = 6 cos φ

Ao = 4.04;

1 3

(

)

φ = cos −1 6W ; W =

[B-45]

VijV jk Vki ~ 1 ⎛⎜ ∂u j ∂u i V V V V ; ; = = + ~ ij ij ij 2 ⎜⎝ ∂xi ∂x j V3

- 179 -

⎞ ⎟ ⎟ ⎠

[B-46]

Appendix B

It can be seen that C µ is a function of the mean strain and rotation rates, the angular velocity of the system rotation, and the turbulence fields.

B-2.3.3. Model constants The model constants have been established to ensure that the model performs well for certain canonical flows. The model constants are:

C 1ε = 1.44;

(Prt )κ

C 2 = 1.9;

= 1.0;

(Prt )ε

= 1.2;

[B-47]

B-2.4. Modeling the turbulent production in the κ − ε models The term Gκ , representing the production of turbulence kinetic energy, is modeled identically for the standard, RNG, and realizable κ − ε models. From the exact equation for the transport of , this term may be defined as

Gκ = − ρ u i′u ′j

∂u j

[B-48]

∂xi

To evaluate Gκ in a manner consistent with the Boussinesq hypothesis,

Gκ = µ tV 2

[B-49]

where V is the modulus of the mean rate-of-strain tensor, defined by Eqs. [B-39] and [B-46].

B-2.5. Standard wall functions Wall functions are a collection of semi-empirical formulas and functions that in effect "bridge'' or "link'' the solution variables at the near-wall cells and the corresponding quantities on the wall. The wall functions comprise: •

laws-of-the-wall for mean velocity and temperature (or other scalars)

•

formulas for near-wall turbulent quantities

The standard wall functions are based on the proposal of Launder and Spalding (1974), and have been most widely used for industrial flows.

B-2.5.1. Momentum The law-of-the-wall for mean velocity yields

U* =

(

1 ln 9.793 y * 0.4187

)

[B-50]

- 180 -

Appendix B

1

U ≡ *

1

U P C µ 4κ P 2

[B-51]

τw / ρ 1

1

ρC µ 4κ P 2 y P y ≡ µ *

[B-52]

The logarithmic law for mean velocity is known to be valid for 30 < y * < 300. In the FV code, the log-law is employed when y * > 11.225. When the mesh is such that y * < 11.225 at the wall-adjacent cells, the FV code applies the laminar stress-strain relationship that can be written as

U * = y*

[B-53]

It should be noted that, in the FV code, the laws-of-the-wall for mean velocity and temperature are based on the wall unit, y * , rather than y + . These quantities are approximately equal in equilibrium turbulent boundary layers. y + is defined as:

y+ =

ρuτ y µ

[B-54]

B-2.5.2. Energy Reynolds' analogy between momentum and energy transport gives a similar logarithmic law for mean temperature. As in the law-of-the-wall for mean velocity, the law-of-the-wall for temperature comprises the following two different laws: •

linear law for the thermal conduction sublayer where conduction is important

•

logarithmic law for the turbulent region where effects of turbulence dominate conduction

The thickness of the thermal conduction layer is, in general, different from the thickness of the (momentum) viscous sublayer, and changes from fluid to fluid. For example, the thickness of the thermal sublayer for a high-Prandtl-number fluid (e.g., oil) is much less than its momentum sublayer thickness. For fluids of low Prandtl numbers (e.g., liquid metal), on the contrary, it is much larger than the momentum sublayer thickness. The law-of-the-wall implemented in the FV code has the following composite form:

- 181 -

Appendix B

(Tw − TP )ρC p C µ 4κ P 2 1

T* =

1

q&

1 1 ⎧ C µ 4κ P 2 2 1 * ⎪Pr y + ρ Pr UP 2 q& ⎪ ⎪⎪ ⎡ 1 ⎤ * = ⎨ Prt ⎢ ln 9.793 y + P ⎥ + κ ⎦ ⎪ ⎣ 1 1 ⎪ C 4κ 2 ⎪ 1 ρ µ P Prt U P2 + (Pr − Prt )U c2 ⎪⎩ 2 q&

(

*

< yT*

)

)

{

where

(y

[B-55]

} (y

*

> yT*

)

P is computed by using the formula given by Jayatilleke (1969):

⎡⎛ Pr P = 9.24⎢⎜⎜ ⎢⎝ Prt ⎣

⎞ ⎟⎟ ⎠

3

4

⎤ − 0.007 Pr ⎡ Prt ⎤ − 1⎥ ⎢1 + 0.28e ⎥⎦ ⎥⎣ ⎦

[B-56]

The non-dimensional thermal sublayer thickness, y T* , in Eq. [B-54] is computed as the y * value at which the linear law and the logarithmic law intersect, given the molecular Prandtl number of the fluid being modeled. The procedure of applying the law-of-the-wall for temperature is as follows. Once the physical properties of the fluid being modeled are specified, its molecular Prandtl number is computed. Then, given the molecular Prandtl number, the thermal sublayer thickness, y T* , is computed from the intersection of the linear and logarithmic profiles, and stored. During the iteration, depending on the y * value at the near-wall cell, either the linear or the logarithmic profile in Eq. [B-55] is applied to compute the wall temperature or heat flux (depending on the type of the thermal boundary conditions).

B-2.5.3. Species When using wall functions for species transport, the FV code assumes that species transport behaves analogously to heat transfer. Similarly to Eq. [B-54], the law-of-the-wall for species can be expressed for constant property flow with no viscous dissipation as

Y ≡ *

(Y

i,w

− Yi )ρC µ 4 κ P 2 1

J i,w

1

⎧Sc y * ⎪ =⎨ ⎡ 1 ⎤ * ⎪Sc t ⎢⎣ 0.4187 ln 9.793 y + Pc ⎥⎦ ⎩

(

)

(y (y

) >y )

*

< y c*

*

* c

[B-57]

where Yi is the local species mass fraction, Sc and Sct are molecular and turbulent Schmidt numbers, and J i , w is the diffusion flux of species i at the wall. Note that Pc and y c* are calculated in a similar way as P and y T* (Eq. [B-56]), with the difference being that the Prandtl numbers are always replaced by the corresponding Schmidt numbers.

- 182 -

Appendix B

B-2.5.4. Turbulence In the κ − ε models the κ equation is solved in the whole domain including the wall-adjacent cells. The boundary condition for κ imposed at the wall is

∂κ =0 ∂n

[B-58]

where n is the local coordinate normal to the wall. The production of kinetic energy, Gκ , and its dissipation rate, ε , at the wall-adjacent cells, which are the source terms in the κ equation, are computed on the basis of the local equilibrium hypothesis. Under this assumption, the production of κ and its dissipation rate are assumed to be equal in the wall-adjacent control volume. Thus, the production of κ is computed from

Gκ ≈ τ w

τw ∂U =τw 1 1 ∂y 0.4187 ρC µ 4 κ P 2 y P

[B-59]

and ε is computed from 3

εP =

3

C µ 4κ P 2

[B-60]

0.4187 y P

The ε equation is not solved at the wall-adjacent cells, but instead is computed using Eq. [B60]. Note that, as shown here, the wall boundary conditions for the solution variables, including mean velocity, temperature, species concentration, κ , and ε , are all taken care of by the wall functions. Therefore, you do not need to be concerned about the boundary conditions at the walls. The standard wall functions work reasonably well for a broad range of wall-bounded flows. However, they tend to become less reliable when the flow situations depart too much from the ideal conditions that are assumed in their derivation. Among others, the constant-shear and local equilibrium hypotheses are the ones that most restrict the universality of the standard wall functions. Accordingly, when the near-wall flows are subjected to severe pressure gradients, and when the flows are in strong non-equilibrium, the quality of the predictions is likely to be compromised.

B-2.6. Convective heat and mass transfer modeling in the κ − ε models Turbulent heat transport is modeled using the concept of Reynolds' analogy to turbulent momentum transfer. The "modeled'' energy equation is thus given by the following:

- 183 -

Appendix B

⎛ ∂ (ρE ) + ∂ [u i (ρE + p )] = ∂ ⎜⎜ k eff ∂T + u i (τ ij )eff ∂t ∂xi ∂x j ⎝ ∂x j

⎞ ⎟ + Sh ⎟ ⎠

[B-61]

( )eff

where k eff is the effective thermal conductivity, E is the total energy, and τ ij

is the

deviatoric stress tensor, defined by Eq. [B-18]

( )eff

The term involving τ ij

represents the viscous heating. Additional terms may appear in the

energy equation, depending on the physical models used. See Section 2.1.3 for more details. For the standard and realizable κ − ε models, the effective thermal conductivity is given by

k eff = k +

C p µt

[B-62]

Prt

where k is the thermal conductivity. The default value of the turbulent Prandtl number is 0.85. For the RNG κ − ε model, the effective thermal conductivity is

k eff = α C p µ eff

[B-63]

where α is calculated from Eq. [B-29], but with α 0 = 1/Pr = k / C p µ . The fact that α varies with µ mol / µ eff , as in Eq. [B-29], is an advantage of the RNG κ − ε model. It is consistent with experimental evidence indicating that the turbulent Prandtl number varies with the molecular Prandtl number and turbulence (Kays, 1994). Eq. [B-29] works well across a very broad range of molecular Prandtl numbers, from liquid metals (Pr ≈ 10-2) to paraffin oils (Pr ≈ 103), which allows heat transfer to be calculated in low-Reynolds-number regions. Eq. [B-29] smoothly predicts the variation of effective Prandtl number from the molecular value ( α = 1/Pr) in the viscosity-dominated region to the fully turbulent value ( α = 1.393) in the fully turbulent regions of the flow. Turbulent mass transfer is treated similarly. For the standard and realizable κ − ε models, the default turbulent Schmidt number is 0.7. For the RNG model, the effective turbulent diffusivity for mass transfer is calculated in a manner that is analogous to the method used for the heat transport. The value of α 0 in Eq. [B-29] is α 0 = 1/Sc, where Sc is the molecular Schmidt number.

B-3. The standard κ − ω model The standard κ − ω model is an empirical model based on model transport equations for the turbulence kinetic energy ( κ ) and the specific dissipation rate ( ω ), which can also be thought of as the ratio of κ to ε (Wilcox, 1998a, 1998b).

- 184 -

Appendix B

B-3.1. Transport equations The turbulence kinetic energy, κ , and the specific dissipation rate, ω , are obtained from the following transport equations:

⎛ ∂ (ρκ ) + ∂ (ρ κ u i ) = ∂ ⎜⎜ Γκ ∂κ ∂t ∂xi ∂x j ⎝ ∂x j

⎞ ⎟ + Gκ − Yκ + S κ ⎟ ⎠

[B-64]

⎛ ∂ (ρω ) + ∂ (ρ ω u i ) = ∂ ⎜⎜ Γω ∂ϖ ∂t ∂xi ∂x j ⎝ ∂x j

⎞ ⎟ + Gω − Yω + S ω ⎟ ⎠

[B-65]

In these equations, Gκ represents the generation of turbulence kinetic energy due to mean velocity gradients. Gϖ represents the generation of ϖ . Γκ and Γω represent the effective diffusivity of κ and ω , respectively. Yκ and Yω represent the dissipation of κ and ω due to turbulence. All of the above terms are calculated as described below. S κ and S ω are userdefined source terms.

B-3.2. Modeling the effective diffusivity The effective diffusivities for the κ − ω model are given by

Γκ = µ +

Γω = µ +

µt

[B-66]

(Prt )κ µt

[B-67]

(Prt )ω

The turbulent viscosity is computed by combining κ and ω as follows:

µt = α *

ρκ ω

[B-68]

B-3.2.1. Low-Reynolds-number correction The coefficient α * damps the turbulent viscosity causing a low-Reynolds-number correction. It is given by

⎛ α * + Re t ⎜ 0 Rκ α * = α ∞* ⎜ ⎜⎜ 1 + Re t Rκ ⎝

⎞ ⎟ ⎟ ⎟⎟ ⎠

[B-69]

- 185 -

Appendix B

Re t =

ρκ ; µω

α 0* =

Rκ = 6;

βi 3

;

β i = 0.072;

[B-70]

Note that, in the high-Reynolds-number form of the κ − ω model, α * = α ∞* = 1 .

B-3.3. Modeling the turbulent production B-3.3.1. Production of κ The term Gκ represents the production of turbulence kinetic energy. It is defined in the same way as for the κ − ε model (see Eq. [B-48]).

B-3.3.2. Production of ω The production of ω is given by

Gω = α

ω Gκ κ

[B-71]

where Gκ is given by Eq. [B-48]. The coefficient α is given by

α α = ∞* α

⎛ α + Re t ⎜ 0 Rω ⎜ ⎜⎜ 1 + Re t Rω ⎝

⎞ ⎟ ⎟ ⎟⎟ ⎠

[B-72]

where Rω = 2.95. α * and Re t are given by Eqs. [B-69] and [B-70] respectively. Note that, in the high-Reynolds-number form of the κ − ω model, α = α ∞ = 1 .

B-3.4. Modeling the turbulence dissipation B-3.4.1. Dissipation of κ The dissipation of κ is given by

Yκ = ρβ * f β * κω

fβ*

⎧1 ⎪ 2 = ⎨1 + 680 X κ ⎪1 + 400 X 2 κ ⎩

[B-73]

Xκ ≤ 0 [B-74]

Xκ > 0

- 186 -

Appendix B

1 ∂κ ∂ω ω 3 ∂x j ∂x j

[B-75]

β * = β i* [1 + ζ * F (M t )]

[B-76]

Xκ =

4 ⎛ ⎞ ⎜ 4 + ⎛⎜ Re t Rβ ⎟⎠ ⎜ 15 ⎝ β i* = β ∞* ⎜ 4 ⎞ ⎜ 1 + ⎛⎜ Re t ⎟ ⎜ Rβ ⎠ ⎝ ⎝

ζ * = 1.5;

Rβ = 8;

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

[B-77]

β ∞* = 0.09;

[B-78]

where Re t and F (M t ) are given by Eqs. [B-70] and [B-83] respectively.

B-3.4.2. Dissipation of ω The dissipation of ω is given by

Yω = ρ β f β ω 2

[B-79]

1 + 70 X ω 1 + 80 X ω

[B-80]

fβ =

Xω =

Ω ij Ω jk Vki

[B-81]

(β ω )

3

* ∞

The strain rate tensor, Vij , is defined by Eq. [B-46]. The rotation rate tensor, Ω ij , is defined by Eq. [B-9]. Also,

⎡

β = β i ⎢1 − ⎣

⎤ β i* * ζ F (M t )⎥ βi ⎦

[B-82]

β i* and F (M t ) are defined by Eqs. [B-77] and [B-83], respectively.

B-3.4.3. Compressibility correction The compressibility function, F (M t ) , is given by

⎧0 F (M t ) = ⎨ 2 2 ⎩ M t − M t0

M t ≤ M t0 M t > M t0

[B-83]

- 187 -

Appendix B

M t2 =

2κ ; a2

M t 0 = 0.25;

a = γRT

[B-84]

Note that, in the high-Reynolds-number form of the κ − ω

model, β i* = β ∞* . In the

incompressible form, β * = β i* .

B-3.5. Model constants α ∞* = 1;

α ∞ = 0.52;

Rκ = 6;

Rω = 2.95;

1 β ∞* = 0.09; β i = 0.072; Rβ = 8; 9 ζ * = 1.5; M t 0 = 0.25; σ κ = 2.0 σ ω = 2.0

α0 = ;

[B-85]

B-3.6. Wall boundary conditions The wall boundary conditions for the κ equation in the κ − ω models are treated in the same way as the equation is treated when enhanced wall treatments are used with the κ − ε models. This means that all boundary conditions for wall-function meshes will correspond to the wall function approach, while for the fine meshes, the appropriate low-Reynolds-number boundary conditions will be applied. The value of ω at the wall is specified as

ρ (u * ) + ωw = ω µ 2

[B-86]

⎛

ω + = min⎜ ω w+ , ⎜ ⎝

⎧⎛ 50 ⎞ 2 ⎪⎜⎜ + ⎟⎟ ⎪ ω w+ = ⎨⎝ κ s ⎠ ⎪ 100 ⎪⎩ κ s+ ⎛

κ s+ = max⎜⎜1.0, ⎝

6

β i (y + )

2

⎞ ⎟ ⎟ ⎠

[B-87]

κ s+ < 25 [B-88]

κ ≥ 25 + s

ρκ s u * ⎞ ⎟ µ ⎟⎠

[B-89]

and k s is the roughness height. In the logarithmic (or turbulent) region, the value of ω + is

ω+ =

1

β ∞*

+ du turb dy +

[B-90]

- 188 -

Appendix B

which leads to the value of ω in the wall cell as

ω=

u*

[B-91]

β ∞* κy

References Choudhury, D., Kim, S.-E., Flannery, W.S., (1993). Calculation of turbulent separated flows using a renormalization group based k-ε turbulence model. American Society of Mechanical Engineers, Fluids Engineering Division (Publication) FED 149, 177 – 187. Jayatilleke, C., (1969). The Influence of Prandtl number and surface roughness on the resistance of the laminar sublayer to momentum and heat transfer. Progress in Heat and Mass Transfer, 1, 193 – 321. Kays, W.M., (1994). Turbulent Prandtl number. Where are we?. Journal of Heat Transfer, 116, 284 – 295. Launder, B.E. and Spalding, D.B., (1972). Lectures in mathematical models of turbulence. Academic Press, London. Launder, B.E. and Spalding, D.B., (1974). The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering, 3, 269 – 289. Reynolds, W.C., (1987). Fundamentals of turbulence for turbulence modeling and simulation. Lecture notes for Von Karman Institute, agard report No. 755. Shih, T.-H., Liou, W.W., Shabbir, A., Yang, Z., Zhu, J., (1995). New k-ε eddy viscosity model for high Reynolds number turbulent flows. Computers and Fluids, 24, 227 – 238. Spalart, P. and Allmaras, S., (1992). A one-equation turbulence model for aerodynamic flows. Technical Report AIAA-92-0439, American Institute of Aeronautics and Astronautics. Wilcox, D.C., (1998a). Reassessment of the scale-determining equation for advanced turbulence models. AIAA Journal, 26, 1299 – 1310. Wilcox, D.C., (1998b). Multiscale model for turbulent flows. AIAA Journal, 26, 1311 – 1320.

- 189 -

APPENDIX C

DIMENSIONLESS NUMBERS

Physical meaning and applied formulas

Harmonic works, « Dimensionless travels »

Appendix C

C-1. Biot number The Biot number is a dimensionless number used in unsteady-state (or transient) heat transfer calculations. It is named after the French physicist Jean Baptiste Biot (17741862), and relates the heat transfer resistance inside and at the surface of a body. The Biot number can be expressed as:

Bi =

hw ⋅ L kr

[C-1]

Figure C-1. Jean Baptiste Biot

C-2. Eckert number The Eckert number is a dimensionless number used in flow calculations. It expresses the relationship between a flow's kinetic energy and enthalpy, and is used to characterize dissipation. It is named for the late professor Ernst R. G. Eckert (1904 – 2004), known for his work with the early development of jet engines and for discovering ways to increase rocket efficiency. The Eckert number can be expressed as:

u2 Ec = C p ∆T

[C-2]

Figure C-2. Ernst R. G. Eckert

C-3. Euler number The Euler number is a dimensionless value common used for analyzing fluid flow dynamics problems where the pressure difference between two points is an important variable. The Euler Number can be interpreted as a measure of the ratio of pressure forces to inertial forces. It is named for Leonhard Euler (1707-1783), mathematician who first explained role of pressure in fluid flow, formulated basic equations of motion and so-called Bernoulli theorem and introduced the concept of cavitation and the principles of centrifugal machinery. The Euler number can be expressed as: Figure C-3. Leonhard Euler

Eu =

p ρ u2

- 193 -

[C-3]

Appendix C

C-4. Froude number The Froude number is a dimensionless parameter measuring of the ratio of the inertia force on an element of fluid to the weight of the fluid element - the inertial force divided by gravitational force. It is named for William Froude (1810 – 1879), engineer, hydrodynamicist and naval architect who was the first to formulate reliable laws for the resistance that water offers to ships and for predicting their stability. The Froude number can be expressed as:

Fr =

u

[C-4]

g Dp

Figure C-4. William Froude

C-5. Grashof number The Grashof number is a dimensionless number in fluid dynamics which approximates the ratio of the buoyancy force to the viscous force acting on a fluid. It is named after the German engineer Franz Grashof (1826 – 1893). The Grashof number can be expressed as:

GrH =

gβ (Ts − T∞ )D p3

ν2

[C-5]

There is an analogous form of the Grashof number used in cases of natural convection mass transfer problems: Figure C-5. Franz Grashof

gD 3p ∆ρ ⎛ ρ ⎞ ⎜ ⎟ GrD = ρ ⎜⎝ µ ⎟⎠

2

[C-6]

C-6. Mach number The Mach number is defined as a ratio of the speed of flow relative to the speed of sound in the flowing fluid. It can be shown that the Mach number is also the ratio of inertial forces (also referred to aerodynamic forces) to elastic forces. The Mach number is named after Austrian physicist and philosopher Ernst Mach (1838 – 1916), an can be expressed as:

Ma =

u

[C-7]

u sound

- 194 -

Figure C-6. Ernst Mach

Appendix C

C-7. Nusselt number The Nusselt number is a dimensionless number which measures the enhancement of heat transfer from a surface which occurs in a 'real' situation, compared to the heat transfer that would be measured if only conduction could occur. Typically it is used to measure the enhancement of heat transfer when convection takes place. It is named after the German engineer Wilhelm Nusselt (1882 – 1957), pioneer in the study of the basic laws of heat transfer, and can be expressed as:

Figure C-7. Wilhelm Nusselt

Nu =

hD p

[C-8]

kf

C-8. Péclet number In physics, the Péclet number is a dimensionless number relating the forced convection of a system to its heat conduction. It is equivalent to the product of the Reynolds number with the Prandtl number. It is named after the French Physicist Jean Claude Eugene Péclet (1793 – 1857) There are various definitions of the Péclet number. The most typical are as follows:

Pe =

D p uC p ρ kf

= Re ⋅ Pr

[C-9]

Figure C-8. Jean Claude Eugene Péclet

C-9. Prandtl number The Prandtl Number is a dimensionless number approximating the ratio of momentum diffusivity and thermal diffusivity. It is named after Ludwig Prandtl (1875 – 1953), German physicist famous for his works in aeronautics. The Prandtl number can be expressed as:

Pr =

Cpµ kf

Figure C-9. Ludwig Prandtl

- 195 -

[C-10]

Appendix C

C-10. Rayleigh number In fluid mechanics, the Rayleigh number for a fluid is a dimensionless number associated with the heat transfer within the fluid. When the Rayleigh number is below the critical value for that fluid, heat transfer is primary in the form of conduction; when it exceeds the critical value, heat transfer is primarily in the form of convection. The Rayleigh number is named after John William Strutt, 3rd Baron of Rayleigh (aka Lord Rayleigh, 1842 – 1919) and is defined as the product of the Grashof number and the Prandtl number.

Ra =

D 3p ρ 2 gβC p ∆T kf µ

= Gr ⋅ Pr

[C-11]

Figure C-10. Lord Rayleigh

C-11. Reynolds number The Reynolds number is the ratio of inertial forces to viscous forces and is used for determining whether a flow will be laminar or turbulent. It is the most important dimensionless number in fluid dynamics and provides a criterion for determining dynamic similitude. It is named after the Irish fluid dynamics engineer Osbourne Reynolds (1842 – 1912), who proposed it in 1883. Typically it is given as follows for flow through a packed bed: Figure C-11. Osbourne Reynolds

Re =

D p uρ

[C-12]

µ

C-12. Schmidt number The Schmidt number is a dimensionless number approximating the ratio of momentum diffusivity and mass diffusivity, and is used to characterize fluid flows in where there are simultaneous momentum and mass diffusion convection processes. It is named after Ernst Schmidt (1892 – 1975), German scientist, pioneer in the field of engineering thermodynamics, especially in heat and mass transfer. The Schmidt number is the mass-transfer equivalent of the Prandtl number, and it can be expressed as:

Sc =

µ

[C-13]

ρD AB

- 196 -

Figure C-12. Ernst Schmidt

Appendix C

C-13. Sherwood number The Sherwood number is a dimensionless number used in mass-transfer operation. It represents the ratio of length scale to the diffusive boundary layer thickness. It is named after Thomas Kilgore Sherwood (1903 – 1976), American chemical engineer whose primary research area was mass transfer and its interaction with flow and with chemical reaction. The Sherwood number is the mass-transfer equivalent of the Nusselt number, and it can be expressed as: Figure C-13. Thomas Kilgore Sherwood

Sh =

kc D p

[C-14]

D AB

C-14. Stanton number The Stanton number is a dimensionless number which measures the ratio of heat transferred into a fluid to the thermal capacity of fluid. It is used to characterize heat transfer in forced convection flows. The Stanton number is named after Sir Thomas Eduard Stanton (1865 – 1931). Stanton's main field of interest was fluid flow and friction, and the related problem of heat transmission.

St =

h Cpρ u

=

Nu Re ⋅ Pr

[C-15]

Figure C-14. Sir Thomas Eduard Stanton

C-15. Strouhal number The Strouhal Number is a dimensionless value useful for analyzing oscillating, unsteady fluid flow dynamics problems. It represents a measure of the ratio of inertial forces due to the unsteadiness of the flow or local acceleration to the inertial forces due to changes in velocity from one point to an other in the flow field.

Figure C-15. Čeněk Strouhal

The Strouhal number, named after the Czech professor of experimental physics Čeněk Strouhal (1850 – 1922), is proportional to the reciprocal of vortex spacing expressed as number of obstacle diameters and is used in momentum transfer in general and Von Karman vortex streets and unsteady state flow calculations in particular. It is defined in the following form:

Sr =

ωD p

[C-16]

u

- 197 -

Appendix C

C-16. Turbulent numbers In the development of the dimensionless analysis there appears variations of some of the aforementioned dimensionless numbers when turbulent flow is considered (i.e., Reynolds, Prandtl and Schmidt turbulent numbers). In these cases, applied formulas to estimate this numbers are as follows:

Re t =

Prt =

Sc t =

D p uρ

[C-17]

µt C p µt

[C-18]

kt

µt

[C-19]

ρD AB

Turbulence intensity used to estimate used to estimate k t from

µt

µt

was set between 1 - 4 %, and Reynolds analogy was

(White, 1991). For first analysis purposes, Prt was assumed as a

constant value within the bed, justified in the fact that majority of experimental results shows a range of variations between 0.75 < Prt < 2 for air and water (Kays, 1994).

References Kays, W.M., (1994). Turbulent Prandtl number. Where are we? Journal of Heat Transfer, 116, 284 – 295. White, F., (1991). Viscous Fluid Flow; McGraw Hill, New York, pp. 482 – 542.

- 198 -

APPENDIX D

FLUID PROPERTIES

ESTIMATION & CFD IMPLEMENTATION

Jodie Coston, « Green background »

Appendix D

D-1. Density D-1.1. Gas @ atmospheric pressure In the case of a gas at atmospheric pressure it was chosen to define the density using the ideal gas law for an incompressible flow, which is implemented by default in the FV code (Fluent Inc., 2005). In this case, the density is computed as:

p op M w

ρ=

[D-1]

RT

where R is the universal gas constant, Mw is the molecular weight of the gas and pop is the operating pressure. In this form, the density depends only on the operating pressure and not on the local relative pressure field.

D-1.2. Fluid @ P ≥ Pc In the case of a supercritical fluid it was chosen to define the density using a cubic equation of state (EOS). In this case the density is computed as:

p op M w

ρ=

[D-2]

zRT

where z is the compressibility factor obtained from solving a cubic EOS, which in general form can be expressed as:

(

)

(

)

z 3 − 1 + B * − uB * z 2 + A* + wB *2 − uB * − uB *2 z − A* B * − wB *2 − wB *3 = 0

[D-3]

where:

A* =

aP R 2T 2

[D-4]

bP RT

[D-5]

and

B* =

Two cubic EOS were used to estimate the density of a fluid in supercritical conditions during the development of this thesis. Redlich-Kwong EOS (Redlich and Kwong, 1949) was used for estimating the density of a pure substance, and Peng-Robinson EOS (Peng and Robinson, 1976) was applied for estimating the density of a mixture. For these equations, a, b, u and w take the values listed in Table D-1.

- 201 -

Appendix D

Equation RedlichKwong

u 1

w

b

a

0

0.08664 RTc Pc

0.42748 R 2Tc2.5 Pc T 0.5

0.07780 RTc Pc

⎛ ⎛ T ⎞ 0.5 ⎞⎤ 0.45724 R 2Tc2 ⎡ ⎢1 + fω ⎜1 − ⎜ ⎟ ⎟⎥ ⎜ ⎜⎝ Tc ⎟⎠ ⎟⎥ Pc ⎢ ⎝ ⎠⎦ ⎣ fω = 0.37464 + 1.54226ω − 0.26992ω 2

2

PengRobinson Table D-1.

2

-1

Constants for the cubic equations of state used

As commercial CFD codes do not offer real-gas equations of state within their default equation system, these equations have to be created and uploaded to the CFD solver if a supercritical fluid has to be modelled. In order to implement the aforementioned EOS within the CFD solver, user-defined functions (UDF) and user-defined equations (UDE) were used. For the UDF creation, an analytical solution of the EOS was compiled using a build-up C++ compiler (incorporated within the CFD software), and included within the equation system of the CFD solver. This strategy allowed the CFD code to compute the density of the fluid as a function of the static pressure and the temperature in every point of the computational grid. An example of a UDF for estimating the density of CO2 using the Redlich-Kwong EOS would be as follows: %UDF EXAMPLE% %PROPERTY TO BE MODIFIED INTO SOLVER: DENSITY% %FLUID: CO2% %Tr: 304.1 K% %Pr: 7380000 Pa% %EOS: REDLICH-KWONG% #include "udf.h" DEFINE_PROPERTY(cell_density, cell, thread) { real ro; real a; real b; real z; real tr; real pr; real a2; real b2; real r; real q; real f; real g; real c; real d; real e; real fi; real z1; real z2; real z3; real t=C_T(cell, thread);

- 202 -

Appendix D

real p=C_P(cell, thread); tr=t/304.1; pr=p/7380000.; a=6456696.46; b=0.02968; a2=0.42747*pr/pow(tr,(5/2)); b2=0.08664*pr/tr; r=a2*b2; q=pow(b2,2.)+b2-a2; f=((-3.*q)-1.)/3.; g=((-27.*r)-(9.*q)-2.)/27.; c=pow((f/3.),3.)+pow((g/2.),2.); if (c>0.0) { d=pow(((-g/2.)+pow(c,0.5)),(1./3.)); e=pow(((-g/2.)-pow(c,0.5)),(1./3.)); z=d+e+(1./3.); } else { fi=a*cos(-6.75*pow(g,2.)/pow(f,3.)); z1=(2*pow((-f/3.),0.5)*cos(fi/3.))+(1/3); z2=(2*pow((-f/3.),0.5)*cos((fi/3.)+2.0944))+(1./3.); z3=(2*pow((-f/3.),0.5)*cos((fi/3.)+4.1888))+(1./3.); if (z1>z2) { if (z1>z3) z=z1; else z=z3; } else { if (z2>z3) z=z2; else z=z3; } } ro=44.*p/(8314.*z*t); }

return ro;

For situations where the static pressure could be considered as constant along the computed model, a simplified density model could be used. This simplified model consisted on expressing the density of the fluid (at constant pressure) as a function of temperature, fitting the results obtained with the EOS to a 6th degree polynomial expression. The obtained polynomial expression (which can be written in general form as shown in Eq. [D-6]), can then be fed into the CFD solver as an UDE.

ρ = aT 6 + bT 5 + cT 4 + dT 3 + eT 2 + fT + g

[D-6]

A sample of the obtained coefficients for Eq. [D-6] for estimating the density of supercritical CO2 can be seen in Table D-2, and a graphical representation of the previously stated polynomial fitting is shown in Figure D-1. In order to quantify the round-up error associated

- 203 -

Appendix D

with the numerical solver, a simulation for a density varying flow of supercritical CO2 in a Venturi tube was set, and obtained numerical values were compared with the prediction of the EOS, resulting in a good agreement between the numerical results and the EOS estimation (see Figure D-2).

P [MPa]

a

b

8

1.13 x 10-9

c

d

e

f

g

-2.65 x 10-6

2.58 x 10-3

-1.33 x 100

3.87 x 102

-5.98 x 104

3.84 x 106

10-7

10-4

10-1

101

104

-6.48 x 105

10

-4.05 x

12

-1.44 x 10-9

3.34 x 10-6

-3.22 x 10-3

1.65 x 100

-4.72 x 102

7.17 x 104

-4.52 x 106

14

-4.78 x 10-10

1.14 x 10-6

-1.14 x 10-3

6.01 x 10-1

-1.78 x 102

2.78 x 104

-1.81 x 106

16

6.61 x 10-11

-1.22 x 10-7

8.58 x 10-5

-2.70 x 10-2

2.75 x 100

3.46 x 102

-6.86 x 104

18

1.79 x 10-10

-3.98 x 10-7

3.65 x 10-4

-1.77 x 10-1

4.77 x 101

-6.78 x 103

3.99 x 105

20

1.43 x 10-10

-3.24 x 10-7

3.04 x 10-4

-1.51 x 10-1

4.18 x 101

-6.12 x 103

3.72 x 105

22

10-11

10-7

10-4

10-2

101

103

2.50 x 105

Table D-2.

8.89 x

10-10

8.77 x

-2.04 x

-7.80 x

1.94 x

3.63 x

-9.76 x

-9.30 x

2.73 x

1.23 x

-4.05 x

Obtained coefficients for Eq. [D-6] for supercritical CO2 850

3

Density [kg/m ]

700 550 400 250 100 310

Figure D-1.

330

350

370 390 Temperature [K]

410

430

8 MPa

10 MPa

12 MPa

14 MPa

16 MPa

18 MPa

20 MPa

22 MPa

450

Density polynomial fitting for supercritical CO2

from CFD

240

ρ

3

250

[kg/m ]

255

245

235 230 230

235

240

ρ Figure D-2.

[kg/m3]

245

250

255

from EOS

Density parity plot for supercritical CO2 (numerical results vs. EOS)

- 204 -

Appendix D

D-1.3. Mixtures In order to study how the density of a mixture was estimated by the CFD software, a 2dimensional cylinder of 5 m length and 0.2 m diameter was created. From one side of the tube CO2 in laminar flow regime was fed to the model. In order to impose a density gradient along the tube, a flux of 0.01 kg/s of toluene vapor was imposed through the pipe walls. Three temperatures (i.e. 298, 310 and 320 K) and four pressures (i.e. 0.1, 6, 8 and 20 MPa) were taken as study cases. For atmospheric pressure situations, incompressible ideal gas law was used to estimate density. For near-critical and supercritical conditions, Peng Robinson equation of state for density was implemented through a C++ subroutine inside the CFD solver. For each of the simulations it was recorded the C7H8 mass fraction together with the mixture density at the central axis of the tube. The numerical results obtained were compared against the estimation through equations of state (Guardo et al., 2005). Low pressure density was determined by ideal gas law for both of the pure components assuming the flow as incompressible. In such way the density is computed by the CFD solver as shown in Eq. [D-7].

ρ=

P

[D-7]

y R ⋅T ⋅ ∑ i Mw

It was expected that for low pressure both methods would give similar results due to the low interaction forces between the molecules. The results obtained at the different temperatures tested (i.e. 298, 310 and 320K) are very similar to the ones shown in Figure D-3. From those results it was checked that the numerical error of the CFD solver when applying the ideal gas law was neglectable.

3

Mixture density [kg/m ]

1.95 1.9 1.85 1.8 1.75 0

0.02

0.04

0.06

0.08

0.1

Mass fraction C7H8 EOS Estimation

Figure D-3.

CFD Estimation

Density estimation for a binary mixture [C7H8/CO2] at 298 K and 0.1 MPa

In order to study high pressure situations, estimation of the mixture density was checked at 6, 8 and 20 MPa. Peng-Robinson EOS was implemented via UDF in order to obtain the compressibility factors. In such way the density is computed by the CFD solver as shown in Eq. [D-8].

- 205 -

Appendix D

ρ=

P

[D-8]

y z ⋅ R ⋅T ⋅ ∑ i Mw

An example of the results obtained for high pressure situations can be seen in Figure D-4. In general, for all situations simulated, it was noticed that despite the fact that the estimation done by the CFD solver is close to the EOS prediction for low solute mass fractions due to the good pure component density estimation, such a simple mixing rule is not able to predict the mixture critical point or phase changes. For all cases analyzed the divergence between CFD estimated results and the EOS prediction increases as the mass fraction of solute increases. According to this, the greatest numerical errors obtained for the CFD simulations were located at the saturation point (equilibrium conditions at the specified pressure and temperature). Maximum error measured when estimating density was of 6.24%, which is an acceptable value.

Mixture density [kg/m3]

300

275

250

225

200 0

0.003

0.006 Mass fraction C7H8

EOS estimation

Figure D-4.

0.009

0.012

CFD Estimation

Density estimation for a binary mixture [C7H8/CO2] at 320 K and 8 MPa

D-2. Viscosity D-2.1. Gas @ atmospheric pressure In the case of a gas at atmospheric pressure it was chosen to define the density using a viscosity power law, which is implemented by default in the FV code (Fluent Inc., 2005). In this case, the viscosity is computed as:

⎛T µ = µ 0 ⎜⎜ ⎝ T0

⎞ ⎟⎟ ⎠

n

[D-9]

where T0 and µ0 are reference values for temperature and viscosity, respectively. Table D-3 shows the values of the coefficients involved in Eq. [D-9] for CO2 and air at moderated pressures and temperatures.

- 206 -

Appendix D

Fluid

µ0 [Pa· s]

T0 [K]

n

CO2

1.37 x 10-5

273.11

0.79

Air

1.716 x 10-5

273.11

0.666

Table D-3.

Viscosity power law coefficients for CO2 and air

D-2.2. Fluid @ P ≥ Pc In the case of a supercritical fluid it was chosen to define the viscosity following the procedure recommended by Lucas (1980), which first calculates the parameters Z1 and Z2 for the reduced temperature of interest:

(

)

Z 1 = 0.807Tr0.618 − 0.357e −0.449Tr + 0.340e −4.058Tr + 0.018 FP0 FQ0 = µ 0ξ 0

[D-10]

0

where FP and FQ are low pressure polarity and quantum correction factors determined as shown in Eqs. [D-12] and [D-13] respectively. To obtain the first of the aforementioned correction factors, a reduced dipole moment is required. Lucas defines this quantity as

σ r = 1.75216 × 10 − 29

σ 2 Pc

[D-11]

Tc2

0

Then, FP values are found as:

0 ≤ σ r < 0.022

FP0 = 1

F = 1 + 30.55(0.292 − z c )

1.72

0 P

FP0 = 1 + 30.55(0.292 − z c )

1.72

0.96 + 0.1(Tr − 0.7 )

0.022 ≤ σ r < 0.075

[D-12]

0.075 ≤ σ r

0

The factor FQ is used only for the quantum gases He, H2 and D2.

[

FQ0 = 1.22Q 0.15 ⎧⎨1 + 0.00385 (Tr − 12 ) ⎩

2

]

1

Mw

Tr − 12 ⎫⎬ ⎭

[D-13] 0

Where Q = 1.38 for He, 0.76 for H2 and 0.52 for D2. For non-quantum gases, FQ = 1. If (1 < Tr < 40) and (0 < PR ≤ 100), then

⎡ ⎤ aPre Z 2 = µ ξ ⎢1 + ⎥ f d −1 ⎣⎢ bPr + (1 + cPr ) ⎦⎥ 0

[D-14]

- 207 -

Appendix D

where µ follows:

a=

0

ξ

is found from Eq. [D-10]. The values of the constants in Eq. [D-14] are defined as

1.245 × 10 −3 5.1726Tr− 0.3286 e Tr

b = a (1.6553Tr − 1.2723) c=

0.4489 3.0578Tr−37.7332 e Tr

d=

1.7368 2.2310Tr− 7.6351 e Tr

[D-15]

e = 1.3088 0.4489

f = 0.9425e −0.1853Tr

After computing Z1 and Z2, we define

Y=

Z2 Z1

[D-16]

and the correction factors FP and FQ,

FP =

FQ =

1 + (FP0 − 1)Y −3 FP0

(

[D-17]

)[

1 + FQ0 − 1 Y −1 − 0.007(ln Y )

4

]

[D-18]

0 Q

F

Finally, the dense gas viscosity is calculated as

µ=

Z 2 FP FQ

[D-19]

ξ

where

⎛ T ξ = 0.176⎜⎜ 3c 4 ⎝ M w Pc

⎞ ⎟⎟ ⎠

1

6

[D-20]

The aforementioned method for viscosity estimation was implemented into the CFD solver via UDF. An example of a UDF for estimating the viscosity of CO2 using the Lucas method would be as follows:

- 208 -

Appendix D

%UDF EXAMPLE% %PROPERTY TO BE MODIFIED INTO SOLVER: VISCOSITY% %FLUID: CO2% %Tr: 304.1 K% %Pr: 7380000 Pa% %METHOD: LUCAS% #include "udf.h" DEFINE_PROPERTY(cell_viscosity, cell, thread) { real mu; real a; real b; real c; real d; real e; real f; real mu_lp; real tr; real pr; real t=C_T(cell, thread); real p=C_P(cell, thread); tr=t/304.1; pr=p/7380000; a=(0.001245/tr)*exp(5.1726*pow(tr,-0.3286)); b=a*((1.6553*tr)-1.2723); c=(0.4489/tr)*exp(3.0578*pow(tr,-37.7332)); d=(1.7368/tr)*exp(2.2310*pow(tr,-7.6351)); e=1.3088; f=0.9425*exp(-0.1853*pow(tr,0.4489)); mu_lp=((0.807*pow(tr,0.618))-(0.357*exp(-0.449*tr))+(0.340*exp(4.058*tr))+0.018)/39194.92; mu=mu_lp*(1+((a*pow(pr,e))/((b*pow(pr,f))+pow((1+(c*pow(pr,d))), -1)))); }

return mu;

For situations where the static pressure could be considered as constant along the computed model, a simplified viscosity model could be used. This simplified model consisted on expressing the viscosity of the fluid (at constant pressure) as a function of temperature, fitting the results obtained with the theoretical method to a 6th degree polynomial expression. The obtained polynomial expression (which can be written in general form as shown in Eq. [D-21]), can then be fed into the CFD solver as an UDE.

µ = aT 6 + bT 5 + cT 4 + dT 3 + eT 2 + fT + g

[D-21]

A sample of the obtained coefficients for Eq. [D-21] for estimating the viscosity of supercritical CO2 can be seen in Table D-4, and a graphical representation of the previously stated polynomial fitting is shown in Figure D-5. In order to quantify the round-up error associated with the numerical solver, a simulation for a viscosity-varying flow of supercritical CO2 in a Venturi tube was set, and obtained numerical values were compared with the prediction of the Lucas method, resulting in a good agreement between the numerical results and the theoretical estimation (see Figure D-6).

- 209 -

Appendix D

P [MPa]

a

8

2.22 x 10-18

-5.27 x 10-15

5.21 x 10-12

-2.74 x 10-9

8.07 x 10-7

-1.26 x 10-4

8.23 x 10-3

10

-1.01 x

10-18

10-15

10-12

10-10

10-7

10-5

-9.52 x 10-4

12

-4.11 x 10-18

9.30 x 10-15

-8.72 x 10-12

4.33 x 10-9

-1.20 x 10-6

1.77 x 10-4

-1.07 x 10-2

14

-4.77 x 10-18

1.09 x 10-14

-1.03 x 10-11

5.17 x 10-9

-1.45 x 10-6

2.16 x 10-4

-1.33 x 10-2

16

-4.18 x 10-18

9.57 x 10-15

-9.08 x 10-12

4.57 x 10-9

-1.28 x 10-6

1.91 x 10-4

-1.18 x 10-2

18

-3.32 x 10-18

7.60 x 10-15

-7.19 x 10-12

3.60 x 10-9

-1.01 x 10-6

1.49 x 10-4

-9.12 x 10-3

20

-2.52 x

10-18

10-15

10-12

10-9

10-7

10-4

-6.41 x 10-3

22

-1.75 x 10-18

6.47 x 10-5

-3.67 x 10-3

Table D-4.

b 2.10 x

5.71 x

c

3.89 x 10-15

-1.78 x

-5.36 x

d

-3.57 x 10-12

7.80 x

2.66 x

e

1.72 x 10-9

-1.84 x

-7.34 x

f

-4.61 x 10-7

2.18 x

1.07 x

g

Obtained coefficients for Eq. [D-21] for supercritical CO2

Viscosity [Pa· s]

9.E-05

6.E-05

2.E-05 310

330

350 8 MPa 16 MPa

Figure D-5.

370 390 Temperature [K] 10 MPa 18 MPa

410

12 MPa 20 MPa

430

450

14 MPa 22 MPa

Viscosity polynomial fitting for supercritical CO2

Viscosity [Pa· s] from CFD

2.16E-05

2.15E-05

2.14E-05

2.13E-05 2.13E-05

2.14E-05

2.15E-05

2.16E-05

Viscosity [Pa· s] from estimation

Figure D-6.

Viscosity parity plot for supercritical CO2 (numerical results vs. Lucas method)

- 210 -

Appendix D

D-2.3. Mixtures To estimate the viscosity of a mixture (at low or high pressure), a similar test to that described in section D-1.3. was developed. In this test, the pure component viscosity was calculated by the aforementioned method of Lucas (1980), imposed to the CFD solver via UDF or UDE. For estimating the viscosity of the mixture, a mass weighted mixing rule was applied:

µ = ∑ Yi ⋅ µ i

[D-22]

i

Results obtained with the CFD solver were compared against those obtained using Lucas Method for mixtures (i.e., see Figure D-7). The results are almost identical for the three temperatures tested in atmospheric pressure situations. From Figure D-7 it can be seen that the calculations through CFD are almost exactly the same as the ones done with the contrast method, proving as done for the density that the CFD solver fits the results expected with the selected methods for comparison.

Viscosity [Pa· s]

1.57E-05

1.53E-05

1.49E-05

1.45E-05 0

0.025

0.05

0.075

0.1

Mass fraction C7H8 Lucas estimation

Figure D-7.

CFD estimation

Viscosity estimation for a binary mixture [C7H8/CO2] at 310 K and 0.1 MPa

Viscosity [Pa· s]

1.E-03

1.E-04

1.E-05

0

0.025

0.05

0.075

0.1

Mass fraction C7H8 CFD Estimation

Figure D-8.

Lucas Estimation

Viscosity estimation for a binary mixture [C7H8/CO2] at 298 K and 60 MPa

- 211 -

Appendix D

An example of the results obtained for the estimation of viscosity in high pressure situations can be seen in Figure D-8. As happened with density, for all cases analyzed the divergence between CFD estimated results and Lucas estimation increases as the mass fraction of solute increases. According to this, the greatest numerical errors obtained for the CFD simulations were located at the saturation point. Measured errors were lower than 10 %, which is an acceptable range (Guardo et al., 2005).

D-3. Thermal conductivity D-3.1. Gas @ atmospheric pressure In the case of a gas at atmospheric pressure it was chosen to define the thermal conductivity as a constant (independent of pressure and temperature). The values used for the thermal conductivities of the simulated gases were taken from experimental data and specific correlations available in the literature (Reid et al., 1987; Yaws, 1999).

D-3.2. Fluid @ P ≥ Pc In the case of a supercritical fluid it was chosen to define the thermal conductivity using excess thermal conductivity correlations (Stiel and Thodos, 1964).

(

)

(k − k )Γz

5 c

= 1.22 × 10 −2 e 0.535 ρ r − 1

(k − k )Γz

5 c

= 1.14 × 10 −2 e 0.67 ρ r − 1.069

(k − k )Γz

5 c

= 2.60 × 10 −3 e1.155 ρ r + 2.016

0

0

0

(

(

) )

ρ r < 0.5

[D-23]

0.5 < ρ r < 2.0

[D-24]

2.0 < ρ r < 2.8

[D-25]

where

⎡ Tc M w3 N 02 ⎤ Γ=⎢ ⎥ 5 4 ⎣ R Pc ⎦

1

6

[D-26]

The aforementioned method for thermal conductivity estimation was implemented into the CFD solver via UDF. An example of a UDF for estimating the thermal conductivity of CO2 using the thermal excess method would be as follows: %UDF EXAMPLE% %PROPERTY TO BE MODIFIED INTO SOLVER: THERMAL CONDUCTIVITY% %FLUID: CO2% %Tr: 304.1 K% %Pr: 7380000 Pa% %METHOD: THERMAL EXCESS% #include "udf.h" DEFINE_PROPERTY(cell_thermal_conductivity, cell, thread) { real ro=C_R(cell, thread);

- 212 -

Appendix D

real real real real

t=C_T(cell, thread); ko; k; ror;

ror=ro/468.58; ko=-0.007215+(0.00008015*t)+(0.000000005477*pow(t,2.)(0.00000000001053*pow(t,3.); if (ror