CFD Simulations in Mechanically Stirred Tank and Flow Field Analysis: Application to the Wastewater (Sugarcane Vinasse) Anaerobic Digestion

Abstract

Anaerobic digestion is a widely used process for waste treatment and energy production. This natural process takes place in a controlled environment, anaerobic digesters. Mixing is one of the main operating parameters. The understanding of the flows during the agitation of the medium is crucial for the optimization of the process yield. In fact, the mass and heat transfers are enhanced by the agitation. However, the complex biochemical reactions can be inhibited with overly vigorous agitation. A detailed and in-depth understanding of the phenomena occurring during agitation requires modeling studies. In this chapter, we propose a general approach, based on computational fluid mechanics (CFD), to analyze the mechanical mixing of an anaerobic reactor. We apply this work to the anaerobic digestion of the sugarcane vinasse, which is a liquid waste generated during the production of alcohol. The single-phase Reynolds-averaged Navier-Stokes (RANS) simulations of mechanical agitation of Newtonian fluids for different rotational speeds are presented. The equations system is closed with the standard k-epsilon turbulence model. The flow field is analyzed with the velocity profiles, the Q and Lambda2 fields, the pressure and the vorticity.

1. Introduction

Anaerobic digestion is a widely used process for organic waste treatment, such as manure, sludge or biowaste. This natural process is based on the digestion of the material in five biological reactions. Each reaction involves a specific group of micro-organisms. First, the complex polymers are hydrolyzed into water-soluble monomers and oligomers by the hydrolytic and fermentative bacteria. Then, the acidogenesis consists in the transformation of the soluble organic molecules into volatile fatty acids, organic acids, hydrogen, carbon dioxide, alcohols and acetate by acidogenic bacteria [1, 2]. This reaction is followed by the acetogenesis, during which the products of hydrolysis and acidogenesis are converted into acetate and carbon dioxide by the action of acetogenic bacteria [1, 2]. Finally, the methanogenesis, the last step of the degradation process, forms methane through two metabolisms. Acetotrophic methanogens transform acetate into methane and carbon dioxide, while hydrogen methanogens combine hydrogen and carbon dioxide to form methane and water [1, 2]. On a laboratory or industrial scale, this natural process takes place in a controlled environment in anaerobic digesters. The main operating parameters are temperature, agitation, organic load and hydraulic retention time.

Anaerobic digestion is based on a set of biological reactions under the action of various groups of micro-organisms. Therefore, the contact between the organic matter and each group of microorganisms must be guaranteed to provide biogas production. This condition is met in the case of adequate digester agitation. The stirring system is used in order to enhance mass and heat transfers. The stirring can be done by recirculation of gas or leachate or mechanical system. The understanding of the flows during the agitation of the medium is crucial for (1) the comprehension of the impact on biochemical reactions and (2) the optimization of the mixing system. According to the choice of agitation mode, the medium can be well-agitated, homogenous or there might be presence of dead zones or isolated turbulent zones. In this context, the flow analysis by numeric approach is a promising method to identify these zones and two mains issues occur. First, to promote the anaerobic digestion process, it is necessary to limit the volume of dead zones without stirring too vigorously. In fact, in the dead zones, the organic matter could be not digested due to the lack of contact between the substrate and each group of micro-organisms. In this sense studies are being conducted to evaluate the volume and the position of dead zones with the aim of reducing their volumes [3, 4, 5]. Second, in the case of highly mixed zones could imply the destruction of the methanogenic centers and inhibit the biochemical reactions. As a consequence, it is evident that anaerobic digestion yields can be improved by the agitation optimization.

The qualitative description of the effect of flow velocity on the biogas yield is illustrated in Figure 1.

Figure 1.

Qualitative description of the effect of flow velocity on the biogas yield.

As a matter of fact, the properties of the medium are impacted by the mixture type. Consequently, the fluid characterization is an important factor. In this work, we propose first to study the mechanical parameter of the sugarcane vinasse, with the consideration of the rheology of the digestion medium. The sugarcane vinasse is a liquid waste generated from the alcohol production. In the case of liquid digestion, defined by a solids content of less than 15%, the medium may both have a Newtonian or non-Newtonian behavior. The behavior of the fluid has a direct consequence on the flows within the digester. Depending on the waste type, the digestion medium will also have a different viscosity. Viscosity is the term to express the resistance of a fluid to flow and motion. It is assumed that the sugarcane vinasse is Newtonian [6].

In order to model the flows within a digester, it is therefore necessary to know both its technical characteristics (geometry, stirring mode) and the fluid properties. All these choices will have an incidence on the mechanical energy and thus on the expended electrical energy to produce methane.

Consequently, the CFD simulations allow to analyze the flows within a digester of a specific waste mixture, considering its rheological properties. In the literature, the CFD models developed are based on the Reynolds-averaged Navier-Stokes (RANS) equations, the large eddy simulation (LES) as well as direct numerical simulation (DNS). DNS method provides the more precise results by solving directly the Navier-Stokes equations but it is also the most time-consuming method. RANS and LES are specific methods developed for turbulence studies. They are both used by many authors. LES method is more time-consuming than RANS method but provides more accurate outcomes. For example, Fan and al. used RANS method for simulating the hydrodynamic behaviors in a stirred tank [7], and Foukrach et al. used RANS method for studying mixing performance in a gas-mixed anaerobic digester [8].

In the current work, we propose a case study on the vinasse. The RANS approach is used, closed with the standard k-epsilon turbulence model. Simulations for different rotational speeds are carried out in order to analyze the flow field in mechanically stirred digester in function of the rotational speed.

2. Theoretical model

2.1 Assumptions

The anaerobic digestion occurs at the mesophilic regime at 37°C. As previously said, it is assumed that the sugarcane vinasse is Newtonian [6] and the fluid is isothermal and incompressible. The rheological properties of the medium are related to the total solid (TS) content of the substrates. The properties of the substrate are detailed in the following part in the Table 1.

Fluid	Density (kg.m^−3)	Dynamic viscosity (Pa.s)	Sources
Vinasse	1044.69	0.0010097	[6]
Water	1000	0.001
Screened and diluted dairy cow manure	1000	0.0009	[9]

Table 1.

Rheological properties of Newtonian fluids.

In the present study, we consider single-phase turbulent flow.

The sugarcane vinasse from the Rivière du Mât distillery is a recalcitrant waste with a chemical oxygen demand (COD) of 86.70 g_O2.L⁻¹, 6.64% TS content, 4.04% volatile solid content and a pH of 4.84 [10]. The biochemical methane potential is 185.59 NL_CH4.kg_COD⁻¹ [10].

In the next section, we present the geometry used in this work and the mesh generation suitable to the mechanical agitation study.

2.2 Geometry and mesh

This study is carried out on an existing geometry in literature. Based on the study carried out by Wu [11], the geometry consists in a 260 mm height tank (liquid height) with a diameter of 152 mm. The digester bottom is conical with a height of 25.88 mm and a diameter of 41 mm. The tank geometry is shown in Figure 2. Wu [11] validated his model from computer automated radioactive particle tracking (CARPT) experiment developed by Hoffmann et al. [9]. Then, Wu [11] also confronted his results with particle image velocimetry (PIV) experiment carried out by Bugay et al. [12]. The geometry of the impeller is shown in Figure 2.

Figure 2.

Cross section of the pilot geometry with the moving zone and the stationary zone (a) and the impeller geometry (b).

The Lightnin A310 impeller, which is a hydrofoil impeller, is used. Hydrofoil impellers were developed for applications where axial flow is important and low shear is desired [13]. The impeller diameter is 62 mm and the axis diameter is 8 mm. The impeller axis is positioned at a height of 50 mm from the bottom of the digester.

The mixing is modeled using the sliding mesh method. The sliding mesh model is a time-dependent solution approach in which the grid surrounding the rotating component(s) physically moves during the solution [13]. Two zones are defined: the stationary zone and the moving zone. The impeller is located in the moving zone. The interface between the two zones is a regular surface mesh size. The interface meshes of the two zones must be identical.

2.3 Governing equations

The sliding mesh (SM) method is suitable for unsteady flow [14]. In our case study, we are in the case of unsteady flow. Indeed, initially, the mechanical agitation is zero and increases gradually until the desired agitation. The flow is steady only when the rotational speed of the stirrer is reached. Therefore, the SM method is used in this work. In order to model the impeller rotation, the digester tank is thus divided in two domains, the rotating zone which contains the impeller and the stationary zone. The arbitrary mesh interface (AMI) is used to link the two domains. The AMI interface is a pair of detached surfaces, giving the AMI1 and AMI2 boundaries. One belongs to the mobile zone and the other to the stationary zone. These two boundaries are identical and have the same initial and boundary conditions. The rotational speed is assigned directly to the moving zone.

The mathematical model is based on RANS equations associated with the forces linked to the impeller rotation: the Coriolis and centrifugal forces. The governing equations for the rotating zone and the stationary zone are respectively:

$$\left\{\begin{array}{c}\frac{\partial u_R}{\partial t}+\nabla \cdot \left(u_R\otimes u_I\right)=-\nabla p+\nabla \cdot \left({\nu }_{\mathit{eff}}\left(\nabla u_I+{\left(\nabla u_I\right)}^T\right)\right)\\ \nabla \cdot u_R=0\end{array}\right.$$E1

$$\left\{\begin{array}{c}\frac{\partial u_I}{\partial t}+\nabla \cdot \left(u_I\otimes u_I\right)=-\nabla p+\nabla \cdot \left({\nu }_{\mathit{eff}}\left(\nabla u_I+{\left(\nabla u_I\right)}^T\right)\right)\\ \nabla \cdot u_I=0\end{array}\right.$$E2

where $u_I$ and $u_R$ are the absolute velocities viewed respectively from the stationary and the rotating frames (m.s⁻¹), $t$ is the time (s), ${\nu }_{\mathit{eff}}$ is the effective kinematic viscosity (m².s⁻¹) and $p$ is the pressure (Pa).

The incompressible solvers on OpenFoam use the modified pressure (kinematic pressure) $P$ (m².s^{− 2}), calculated by the following equation:

$$P=\frac{p}{\rho }$$E3

where $p$ is the pressure (kg.m⁻¹.s⁻²) in and $\rho$ the fluid density (kg.m⁻³).

The system is closed with the standard $k-\epsilon$ turbulence model, where $k$ is the turbulent kinetic energy (m².s⁻²) and $\epsilon$ is the turbulent dissipation (m².s⁻³). The model implemented in OpenFOAM does not include the buoyancy contribution. The two equations of this model are:

$$\frac{\partial }{\partial t}\left(k\right)=\nabla \cdot \left(D_k\nabla k\right)+G_k-\frac{2}{3}\left(\nabla \cdot u\right)k-\epsilon +S_k$$E4

$$\frac{\partial }{\partial t}\left(\epsilon \right)=\nabla \cdot \left(D_{\epsilon }{\nabla }_{\epsilon }\right)+\frac{C_1{G}_k\epsilon }{k}-\left(\frac{2}{3}{C}_1-C_{3,\mathit{RDT}}\right)\left(\nabla \cdot u\right)\epsilon -C_2\frac{{\epsilon }^2}{k}+S_{\epsilon }$$E5

where $u$ is equal to $u_I$ or $u_R$, $S_k$ is the internal source term for $k$, $S_{\epsilon }$ is the internal source term for $\epsilon$, $G_k$ is the production of $k$, $C_1$ and $C_2$ are model constants, $D_k$ the effective diffusivity for $k$ and $D_{\epsilon }$ is the effective diffusivity for $\epsilon$. The values of the model constants used in OpenFoam, as in the standard model are: $C_{\mu }=0.09,$$C_1=1.44,$$C_2=1.92,$${\sigma }_k=1$ and ${\sigma }_{\epsilon }=1.3$ [15].

The production term $G_k$ is:

$$G_k={\nu }_t{S}^2$$E6

The assumption of the standard $k-\epsilon$ turbulence model turbulence is the following expression of the turbulent viscosity ${\nu }_t$ (m².s⁻¹):

$${\nu }_t=C_{\mu }\frac{k^2}{\epsilon }$$E7

For the initialization of the numerical simulations, the turbulent kinetic energy for isotropic turbulence $k_0$ and the turbulent dissipation rate ${\epsilon }_0$ are calculated as follows:

$$k_0=\frac{3}{2}{\left(I\left|u_{\mathit{ref}}\right|\right)}^2$$E8

$${\epsilon }_0=\frac{{C_{\mu }}^{0.75}{k}_0^{1.5}}{L}$$E9

where $I$ (%) is the turbulence intensity with a default value of 0.05, $u_{\mathit{ref}}$ (m.s⁻¹) is a reference velocity (the stirring velocity) and $L$ (m) is the reference length scale (the digester radius).

The mesh must be refined at the walls. In fact, the standard turbulent model is derived under the assumption of a high local turbulent Reynolds number [16]. This low turbulent Reynolds number region is called the viscous sub layer. Launder and Spalding suggested the following wall function equation [15] to reduce cell number in these zones:

$$u^{+}=\frac{1}{\kappa }ln{y}^{+}+C$$E10

where $u^{+}=\left|u\right|/u_T$ is the dimensionless velocity where $u_T$ is the friction velocity (or shear velocity). $y^{+}=u_{T}y/\nu$ is the dimensionless local Reynolds number where $y$ is the width of the boundary layer (normal distance from the wall). The range of the local Reynolds number is $11.06\le y^{+}\le 300$ [17]. Platteeuw et al. recommend that the local Reynolds number of the first cell should be in the range of $20\le y^{+}\le 100$ [16]. $\kappa$ is the von Karman constant and $C$ is a parameter related to the wall roughness. In the case of smooth walls, the values are respectively 0.419 and 5.24.

The boundary conditions are:

$$\begin{array}{ccc}\mathrm{No} \text{slip}:u=0,& \frac{\partial }{\partial n}p=0,& \mathrm{on} {\mathrm{\Gamma }}_0\\ \text{Imposed velocity field}:v_r=\mathit{\omega R},& \frac{\partial }{\partial n}p=0,& \mathrm{on} {\mathrm{\Gamma }}_1\end{array}$$E11

where the index $r$ is the center of the elementary surfaces of each cell the agitator, $\omega$ is the angular velocity, $R$ is the impeller radius, ${\mathrm{\Gamma }}_0$ is the digester walls and ${\mathrm{\Gamma }}_1$ is the impeller walls.

The rotational speed of the agitator is defined by the speed of rotation of the mobile zone of the mesh. A sudden increase in the stirring speed at the first time step causes a sharp increase in the CFL. Thus, a velocity ramp is necessary for the simulations starting to maintain the Courant number lower than 0.5.

2.4 Rheological expression

We refer to the pseudo-plastic model of Ostwald. It allows to express the shear stress $\tau$ (Pa) as a function of the shear rate $\dot{\gamma }$ (s⁻¹). It is calibrated by two parameters: the consistency index $k$ (Pa.sⁿ) and the flow index $n$. It is expressed as follows:

$$\tau =k{\dot{\gamma }}^n$$E12

The viscosity is expressed as follows:

$$\eta =k{\dot{\gamma }}^{n-1}$$E13

The components of the tensor of the shear rate are:

$${\dot{\gamma }}_{\mathit{ij}}=\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$$E14

$$S_{\mathit{ij}}=\eta \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)=k{\dot{\gamma }}^{n-1}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$$E15

The rheological properties of the vinasse, the water and diluted dairy cow manure are shown in Table 1. It can be seen that the properties of the vinasse are closed to the properties of water. The difference between these two substrates is mainly the presence of suspensions. In anaerobic digestion process, the particles settle, accumulate at the bottom of the digester and result in the formation of sludge. In the present study, as a monophasic model is considered, the suspensions are not taken in account.

2.5 Flow field analysis

The Reynolds number for mechanical agitation of Newtonian fluids is [13, 18, 19]:

$$Re=\frac{\rho {\mathit{Nd}}^2}{\mu }$$E16

where $N$ is the impeller speed (rev.s⁻¹ or rps), $\rho$ is the fluid density (kg.m⁻³), d is the impeller diameter (m) and $\mu$ is the dynamic viscosity (Pa.s).

The generalized Reynolds number valid for non-Newtonian fluids is [5]:

$${Re}_g=\frac{\rho U_{\infty }^{2-n}{d}^n}{K{\left(0.75+\frac{0.25}{n}\right)}^n{8}^{n-1}}$$E17

where $U_{\infty }$ is the average velocity of the fluid (m.s⁻¹), $k$ is the consistency index (Pa.sⁿ), $n$ is the flow index and $d$ is the reference length (m).

The evaluation of the energy consumption and consequently its economic impact is possible with the calculation of the power calculation. The power consumption $P$ (W) calculation based on the torque of the impeller is:

$$P=2\mathit{\pi NT}$$E18

where $N$ is the impeller rotational speed and $T$ is the impeller torque (N.m). The torque is defined as [20]:

$$T=\left(\underset{S}{\int }\left(\overrightarrow{r}\ast \left(\overline{\overline{\tau }}\cdot \widehat{n}\right)\right)\mathit{dS}\right)\cdot \widehat{\alpha }$$E19

where $S$ represents the surfaces including the rotating parts, $\overline{\overline{\tau }}$ Is the total stress tensor, $\widehat{n}$ is a unit vector normal to the surface, $\overrightarrow{r}$ is the position vector and $\widehat{\alpha }$ is a unit vector parallel to the rotation axis.

The mixing energy level (MEL) (W.m⁻³) is calculated by [21]:

$$\mathit{MEL}=\frac{P}{V}$$E20

where $P$ is the power consumption (W) and $V$ is the working volume (m³).

The flow rate through the moving zone $Q$ (m³.s⁻¹) is used for the evaluation of the fluid circulation through the digester. We create a surface for the discharge zone for computing the flow rate. The surface consists of two discs above and below the stirrer for the axial flow and a cylinder (AMI section) for the radial flow.

The power and flow (also called pumping number) numbers are two dimensionless numbers commonly used to characterize the stirred tank flows and mixing processes. The flow number is a measure of the pumping capacity of an impeller [13]. The power number is a dimensionless parameter that provides a measure of the power requirements for the operation of an impeller [13]. These two mixing parameters can be calculated from CFD results and allow to compare the results of simulations with the experimental results. In our case, the experimental values for hydrofoil impeller are $N_P$ equal to 0.3 and $N_Q$ equal to 0.56 [11, 22]. The power and flow numbers are respectively expressed by:

$$N_P=\frac{P}{\rho N^3{d}^5}$$E21

$$N_Q=\frac{Q}{{\mathit{Nd}}^3}$$E22

The spatially average characteristic velocity gradient $G$ (s⁻¹) is used for describing the rotational speed through the digester [23], as well as characterizing the turbulent shear rate [22]. It is defined with the following expression [22, 23]:

$$G={\left(\frac{\overline{\epsilon }}{\nu }\right)}^{1/2}$$E23

where $\nu$ is the kinematic viscosity and $\overline{\epsilon }$ (m².s⁻³) is the global average turbulent energy dissipation rate calculated as follow [22, 23]:

$$\overline{\epsilon }=\left(\frac{N_P{N}^3{d}^5}{V}\right)$$E24

The circulation time $t_c$ is expressed by [22, 24]:

$$t_c=\frac{V}{N_Q{\mathit{Nd}}^3}=\frac{V}{Q}$$E25

The vorticity $\xi$ (s⁻¹) is a measure of the local rotation in the fluid. More specifically, it is a vector field that gives a microscopic measure of the rotation at any point in the fluid. The vorticity describes the local spinning motion. The helicity $H$ provides indications on the alignment of the vorticity vector and the velocity vector $U$. It allows to illustrate the longitudinal vortices, or spiral motion, as is often found in vortex cores.

The vorticity is defined as the curl of the velocity vector:

$$\xi =\nabla \times U$$E26

The helicity is expressed as:

$$H=U\cdot \xi =U\cdot \left(\nabla \times U\right)$$E27

The angle between the vorticity vector and the velocity vector (which is 0° or 180° in a vortex core) is given by:

$$\alpha ={cos}^{-1}\left(\frac{H}{\left|\xi \right|\left|U\right|}\right)$$E28

2.6 CFD simulations

The numerical simulations are performed on OpenFOAM software. The PIMPLE algorithm is used, which is a combination of the PISO (pressure-implicit with splitting of operators) [25, 26] and SIMPLE (semi-implicit method for the pressure linked equations) [27] algorithms. The PISO algorithm is pressure implicit with splitting of operators [25] and allows the use of the sliding mesh method. The SIMPLE algorithm is the semi-implicit method for pressure linked equations [28]. The pimpleFoam solver assumes incompressible, unsteady and viscous flows. The simulations are launched in parallel on 16 processors. The simulation time, for which we get a steady state, is a few seconds. The calculation time is about a few days.

3. Results and discussions

3.1 Mesh

The tank geometry and the mesh are done on OpenFOAM and the impeller geometry on Sketch Up. Then, the OpenFOAM snappyHexMesh command is used to obtain the final mesh. The mesh is refined to meet the quality criteria of mesh such as cell skewness and mesh non-orthogonality. To meet these criteria, we reach a very large number of cells and thus long computing times. Consequently, we do not perform grid independence tests.

The surface mesh of the digester tank and the cross section of the mesh at the impeller level are shown in Figures 3 and 4.

Figure 3.

Surface mesh of the digester: top (left) and wall (right).

Figure 4.

Cross section of the mesh at the impeller level.

The characteristics of the mesh are shown in Table 2. In total, this mesh is composed of 2,841,535 cells with mostly hexahedra. The total mesh faces is 9,150,640. The impeller edge has the larger amount of faces with a value of 465,732 due to its geometry. Concerning the interface between the moving zone and the stationary zone (AMI1 and AMI2), the two surface meshes must be identical with the same faces number in order to avoid numerical errors. Globally, the maximum cell skewness is 3.4 and the maximum mesh non-orthogonality is 65.0, which reflect an acceptable mesh quality. The dimensionless local Reynolds number $y^{+}$ is below 1.

Characteristics
Total mesh cells	2,841,535
Total mesh points	3,485,073
Total mesh faces	9,150,640
Total mesh internal faces	8,591,189
Number of hexahedra	2,578,249
Number of prisms	10,298
Number of wedges	3
Number of polyhedra	252,985
Maximum cell skewness	3.40025
Maximum mesh non-orthogonality	65.0182
Average mesh non-orthogonality	21.3783
Breakdown of polyhedra by number of faces
Faces	Number of cells
<10	216,780
≥10	36,205
Patch topology
Patch	Faces	Points
Upper edge (moving zone)	5719	5952
Bottom edge (moving zone)	4800	4881
Impeller edge (moving zone)	465,732	483,277
Upper edge (stationary zone)	3200	3360
Bottom edge (stationary zone)	3200	3360
Lateral wall (stationary zone)	25,760	25,760
AMI1	25,760	25,760
AMI2	25,760	25,760

Table 2.

Mesh characteristics.

3.2 Velocity profiles and shear stress

In this part, we present the simulation results for four rotational speeds: 20, 40, 60 and 100 rpm. Figure 5 shows the cross section of the velocity field within the digester at the agitator height and Figure 6 extends the longitudinal section of the velocity flied.

Figure 5.

Axial cross-section of the velocity flied (m.s⁻¹) within the digester at the agitator level (z = 0.06 m) for four rotational speeds: 20, 40, 60 and 100 rpm (the maximum velocity in m.s⁻¹ at the impeller extremity is mentioned below the figures).

Figure 6.

Longitudinal cross-section of the velocity field (m.s⁻¹) within the digester for four rotational speeds: 20, 40, 60 and 100 rpm (the maximum velocity in m.s⁻¹ at the impeller extremity is mentioned below the figures).

The results obtained are consistent. In fact, by increasing the rotation velocity of the stirrer, the maximum velocity increases, the volume of dead zones decreases and the volume of continuously agitated zones increases. The maximum velocity value obtained for each case is the agitator velocity. The maximum velocity is obtained at the blade extremities. The velocity decreases gradually with the distance from the agitator.

For each case, we find that the velocity at the edges of the digester is zero. Likewise, the velocity at the axis of the agitator increases with the increase of the rotational speed. The interest of the study of the velocity field at this level relates in particular to the influence of the agitation on the biofilm which tends to develop on the axis of the agitator. The interest is even greater in the case of intermittent mixing. This mode of agitation is particularly interesting because it leads to the reduction of energy consumption. Overall, the optimization of the mechanical agitation must allow to properly stir the medium, so as to improve the biogas yields, while reducing the rotational speed and frequency of agitation.

From both the axial and longitudinal sections of the flow field, we observe that the maximum velocity field is mainly concentrated in the area close to the agitator and quickly becomes zero away from the influence zone of the impeller blades.

We assume that the dead zones are defined by a velocity magnitude less than 0.001 m.s⁻¹ as [5]. Considering this hypothesis, we find that, in our case study, almost the whole area is dead zones except in the zone of influence of the agitator, which is more important with the increase of the rotation speed. Regarding the contours of the longitudinal cross-section, it can be seen that low-speed velocities are observed near the impeller axis. Furthermore, the area at the upper edges of the digester as well as the lower area of the conical section remains dead zones despite a rotational speed of 100 rpm. CFD simulations are useful for defining the volume of dead zones. Indeed, the dead zones within a digester lead to the heterogeneity of the medium. There will be areas where the substrate will not be digested or acidic areas. This can eventually lead to intoxication of the digester generating significant costs on an industrial scale.

Experimentally it would result in an accumulation of solid and undigested matters, and thus the formation of sludge. It is therefore necessary to adapt the agitation according to the digester technology (with or without sludge bed). In addition, in the case of anaerobic digestion of a recalcitrant waste, such as vinasse, it is important to apply sufficient agitation to re-suspend the undigested matter in order to improve yields.

Similarly, since the upper part of the digester is not agitated, the digestion medium would not be homogenized. Overall, this mode of agitation (mechanical, a three-bladed agitator) would result in a significant stratification of the reaction medium with the presence of dead zones.

In view of the obtained results, it would be recommended to use a stirring system with several rows of blades to agitate the lower and upper areas of the digester to avoid stratification of the digestion medium. In addition, it would be desirable that the diameter of the stirrer be close to that of the digester in order to limit the development of biofilm on the walls. In addition, theses simulations provide information on the velocity field within the digester as a function of agitator rotational speed; however, we have no information on their impact on the biochemical reactions and the biogas production and quality.

A study has been carried out on the impact of shear stress and impeller design on the production of biogas in anaerobic digesters [29]. An impact of the shear stress on the biogas production was highlighted [29] and the abrasion of the anaerobic sludge granule due to the shear rate increase above 5 s⁻¹ [29, 30]. Therefore, it was suggested to estimate the shear stress value for the impeller type and mixing rotation [29]. In the present study, the shear stress is negligible in the digester expected at the wall and impeller levels. We obtained 0.31 Pa at 20 rpm and 1.24 Pa at 100 rpm at the tank wall. At the impeller wall, we reported 453 Pa at 20 rpm and 1883 Pa at 100 rpm. Therefore, the shear stress was multiplied by 4.15 at the impeller wall and by 3.9 at the digester wall, with the increase of the stirring speed from 20 to 100 rpm.

CFD simulations accompanied by experimental studies allow to relate the stirring velocity, the volume of dead zones, the homogenization of the medium (physicochemical properties) and the biogas yields. Thus, a better understanding of the physical phenomena involved in the digesters is necessary for the optimization of anaerobic digestion on an industrial scale. An experimental study on the anaerobic digestion of vinasse on this same geometry (digester and agitation system) would provide knowledge on the yields obtained with the velocity fields obtained in modeling. This would make it possible to relate the yields obtained experimentally to the flows.

3.3 Pressure

During the anaerobic digestion process, the pressure applied on the medium has an impact. Indeed, there are processes involving free microorganisms that are in the form of flocs or biofilm [31]. There are two types of biofilms, those formed on a mineral or organic support (fixed or mobile) and granules (natural agglomeration of microorganisms from a few tens of microns to several millimeters in diameter). In fact, the flow study through the digester provide information on the pressure imposed on the digestion medium. The interest is to evaluate the impact of pressure on the flocs and biofilms.

Figures 7 and 8 show the kinematic pressure in continuous flow regime at 40 and 100 rpm respectively. The pressure is due to the impeller rotation. At 20 rpm, it varies from 0 to 5.86.10⁻⁵ m².s⁻² (0 to 5.86.10⁻² Pa). The pressure range at 40 rpm is −0.50 to 0.40 m².s⁻² which corresponds to −500 to 400 Pa. At 60 rpm, the pressure varies from −0.73 to 0.61 m².s⁻² (−730 to 610 Pa). At 100 rpm, it varies from −1.16 to 1.05 m².s⁻² (−1160 to 1050 Pa).

Figure 7.

Kinematic pressure (m².s⁻²) at 40 rpm (z = 0.06 m, y = 0).

Figure 8.

Kinematic pressure (m².s⁻²) at 100 rpm (z = 0.06 m, y = 0).

3.4 Turbulence phenomena analysis

The turbulence phenomena can be studied from different parameters. In the present study, we propose the description of the Reynolds number, the isotropic turbulence, the Q field, the Lambda2 fields and the vorticity.

Considering the impeller diameter for the characteristic length, the Reynolds number is 1326 at 20 rpm, 2651 at 40 rpm, 3977 at 60 rpm and 6629 at 100 rpm. Therefore, local turbulence phenomena can occur around the impeller, however, the turbulence phenomena dissipate rapidly away from the impeller. If the isotropic turbulence k is very low, we consider that we are in a laminar regime. The maximum k value is 2.2 × 10⁻¹⁶ m².s⁻² at 20 and 40 rpm, 3.0 × 10⁻¹¹ m².s⁻² at 60 rpm and 3.9 × 10⁻⁸ m².s⁻² at 100 rpm. Therefore, we observe these local important values at the impeller level.

Q and Lambda2 fields provide a precise description of the turbulence and local rotation. The Lambda2 (s⁻²) function object computes the second largest eigenvalue of the sum of the square of the symmetrical and anti-symmetrical parts of the velocity gradient tensor. Q iso-surfaces are good indicators of turbulent flow structures. The Q function object computes the second invariant of the velocity gradient tensor (s⁻²):

$$Q=\frac{1}{2}\left[{\left(\mathit{tr}\left(\nabla u\right)\right)}^2-\mathit{tr}\left(\nabla u\cdot \nabla u\right)\right]$$E29

Figure 9 presents the velocity flood with the Lambda2 contours. The Lambda2 and velocity profiles are similar in zones between the blades but different at the impeller extremities.

Figure 9.

Velocity (m.s⁻¹) flood and Lambda2 (s⁻²) contour at 40 rpm (z = 0.06 m).

Figure 10 shows the Q field at the impeller level at 40 and 100 rpm. The extremum values of the Q field are obtained at the impeller extremities. The Q field varies from – 5171 to 39,405 s⁻² at 40 rpm and from −27,000 to 289,000 s⁻² at 100 rpm. The value is close to zero in areas away from the impeller. The Q field varies from −1383 to 8924 s⁻² at 20 rpm and from −9314 to 94,137 s⁻² at 60 rpm.

Figure 10.

Q field (s⁻²) at the impeller level at 40 and 100 rpm.

The Lambda2 field at the impeller level at 40 rpm is shown in Figure 11. It varies from – 3677 to 31,177 s⁻² at 40 rpm and from – 15,561 to 226,008 s^{− 2} at 100 rpm. Both the Q and Lambda2 field maximum values are about seven times higher at 100 rpm than at 40 rpm. The Lambda2 field varies from −1125 to 8449 s⁻² at 20 rpm and from −5234 to 73,482 s⁻² at 60 rpm.

Figure 11.

Lambda2 field (s⁻²) at the impeller level at 40 rpm.

Figure 12 presents the vorticity module (s⁻¹) at different digester heights in continuous flow regime at 40 rpm (z = 0.06 m, z = 0.08 m, z = 0.259 m and y = 0). The vorticity is a pseudo vector field that describes the local spinning motion (the curl of the velocity). The maximum vorticity module is 314.8 s⁻¹ at 20 rpm, 621.7 s⁻¹ at 40 rpm, 959.4 s⁻¹ at 60 rpm and 1680.7 s⁻¹ at 100 rpm.

Figure 12.

Vorticity module (s⁻¹) in continuous flow regime at 40 rpm (z = 0.06 m, z = 0.08 m, z = 0.259 m and y = 0).

The extremum values are observed at the impeller level. We notice that the maximum value of vorticity is multiplied by 2.7 from 40 to 100 rpm. We can consequently note that the variations of the Q and Lambda2 fields are more significant than the variation of the vorticity with the increase of the impeller mixing speed.

Figure 13 summarizes the variation of maximum vorticity, pressure, Q and Lambda2 values in function of mixing speed. The Q maximum value variation is the most pronounced among the parameters observed.

Figure 13.

Variation of maximum vorticity, pressure, Q and Lambda2 values in function of mixing speed.

4. Conclusion

Three-dimensional CFD simulations of the mechanically stirred digester were carried out. The model was based on RANS equations and standard k-epsilon turbulence model and the fluid was Newtonian. The sliding mesh method was conducted to characterize the impeller rotation. The outcomes were the velocity profiles, the shear stress, the pressure and the turbulence phenomena analysis with the Q and Lambda2 fields, the Reynolds number, the isotropic intensity and the vorticity. The area of influence of the agitator was in the zone of the blades. The volume of dead zones was important for the four rotational velocities studied. The maximum velocity was observed at the blades extremities. In conclusion, several rows of paddles would be required to reduce the dead zone volume.

In perspective of this work, the flow study will be carried out for different tank geometry, agitation configuration and fluid rheology (non-Newtonian).

Acknowledgments

This work was supported by the Region Reunion (France) as part of the funding of a research thesis in the PIMENT laboratory at the Reunion Island University.