**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 wellagitated, 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–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**.

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].

**Figure 1.** *Qualitative description of the effect of flow velocity on the biogas yield.*

*CFD Simulations in Mechanically Stirred Tank and Flow Field Analysis: Application… DOI: http://dx.doi.org/10.5772/intechopen.93926*

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**.

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 gO2.L<sup>1</sup> , 6.64% TS content, 4.04% volatile solid content and a pH of 4.84 [10]. The biochemical methane potential is 185.59 NLCH4.kgCOD<sup>1</sup> [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


**Table 1.** *Rheological properties of Newtonian fluids.* 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**.

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

#### **Figure 2.**

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

*CFD Simulations in Mechanically Stirred Tank and Flow Field Analysis: Application… DOI: http://dx.doi.org/10.5772/intechopen.93926*

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:

$$\begin{cases} \frac{\partial \boldsymbol{u}\_{R}}{\partial t} + \nabla \cdot (\boldsymbol{u}\_{R} \otimes \boldsymbol{u}\_{I}) = -\nabla p + \nabla \cdot \left(\nu\_{\text{eff}} \left(\nabla \boldsymbol{u}\_{I} + \left(\nabla \boldsymbol{u}\_{I}\right)^{T}\right)\right) \\ \qquad \nabla \cdot \boldsymbol{u}\_{R} = \mathbf{0} \end{cases} \tag{1}$$

$$\begin{cases} \frac{\partial u\_I}{\partial t} + \nabla \cdot (u\_I \otimes u\_I) = -\nabla p + \nabla \cdot \left(\nu\_{\text{eff}} \left(\nabla u\_I + \left(\nabla u\_I\right)^T\right)\right) \\ \qquad \nabla \cdot u\_I = 0 \end{cases} \tag{2}$$

where *uI* and *uR* are the absolute velocities viewed respectively from the stationary and the rotating frames (m.s�<sup>1</sup> ), *t* is the time (s), *νeff* is the effective kinematic viscosity (m<sup>2</sup> .s�<sup>1</sup> ) and *p* is the pressure (Pa).

The incompressible solvers on OpenFoam use the modified pressure (kinematic pressure) *P* (m<sup>2</sup> .s� <sup>2</sup> ), calculated by the following equation:

$$P = \frac{p}{\rho} \tag{3}$$

where *p* is the pressure (kg.m�<sup>1</sup> .s�<sup>2</sup> ) in and *ρ* the fluid density (kg.m�<sup>3</sup> ).

The system is closed with the standard *k* � *ε* turbulence model, where *k* is the turbulent kinetic energy (m<sup>2</sup> .s�<sup>2</sup> ) and *ε* is the turbulent dissipation (m<sup>2</sup> .s�<sup>3</sup> ). The model implemented in OpenFOAM does not include the buoyancy contribution. The two equations of this model are:

$$\frac{\partial}{\partial t}(k) = \nabla \cdot (D\_k \nabla k) + \mathbf{G}\_k - \frac{2}{3}(\nabla \cdot u)k - \varepsilon + \mathbf{S}\_k \tag{4}$$

$$\frac{d}{dt}(\varepsilon) = \nabla \cdot (D\_\varepsilon \nabla\_\varepsilon) + \frac{C\_1 G\_k \varepsilon}{k} - \left(\frac{2}{3} C\_1 - C\_{3,RT}\right) (\nabla \cdot u) \varepsilon - C\_2 \frac{\varepsilon^2}{k} + S\_\varepsilon \tag{5}$$

where *u* is equal to *uI* or *uR*, *Sk* is the internal source term for *k*, *S<sup>ε</sup>* is the internal source term for *ε*, *Gk* is the production of *k*, *C*<sup>1</sup> and *C*<sup>2</sup> are model constants, *Dk* the effective diffusivity for *k* and *D<sup>ε</sup>* is the effective diffusivity for *ε*. The values of the model constants used in OpenFoam, as in the standard model are: *C<sup>μ</sup>* ¼ 0*:*09, *C*<sup>1</sup> ¼ 1*:*44, *C*<sup>2</sup> ¼ 1*:*92, *σ<sup>k</sup>* ¼ 1 and *σε* ¼ 1*:*3 [15].

The production term *Gk* is:

$$\mathbf{G}\_k = \nu\_t \mathbf{S}^2 \tag{6}$$

The assumption of the standard *k* � *ε* turbulence model turbulence is the following expression of the turbulent viscosity *ν<sup>t</sup>* (m<sup>2</sup> .s�<sup>1</sup> ):

$$
\omega\_t = \mathbf{C}\_{\mu} \frac{k^2}{\varepsilon} \tag{7}
$$

For the initialization of the numerical simulations, the turbulent kinetic energy for isotropic turbulence *k*<sup>0</sup> and the turbulent dissipation rate *ε*<sup>0</sup> are calculated as follows:

$$k\_0 = \frac{3}{2} \left( I|u\_{\text{ref}}| \right)^2 \tag{8}$$

$$
\varepsilon\_0 = \frac{\mathbf{C}\_{\mu}^{0.75} k\_0^{1.5}}{L} \tag{9}
$$

where *I* (%) is the turbulence intensity with a default value of 0.05, *uref* (m.s�<sup>1</sup> ) 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 localturbulent 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:

$$
\mu^{+} = \frac{1}{\kappa} \ln y^{+} + \mathcal{C} \tag{10}
$$

where *u*<sup>þ</sup> ¼ j j *u =uT* is the dimensionless velocity where *uT* is the friction velocity (or shear velocity). *y*<sup>þ</sup> ¼ *uTy=ν* 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≤*y*<sup>þ</sup> ≤300 [17]. Platteeuw et al. recommend that the local Reynolds number of the first cell should be in the range of 20≤*y*<sup>þ</sup> ≤100 [16]. *κ* 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{aligned} \text{No slip: } u &= 0, & \frac{\partial}{\partial n} p &= 0, \quad \text{on } \Gamma\_0 \\ \text{Impposed velocity field: } v\_r &= oR, & \frac{\partial}{\partial n} p &= 0, \quad \text{on } \Gamma\_1 \end{aligned} \tag{11}$$

where the index *r* is the center of the elementary surfaces of each cell the agitator, *ω* is the angular velocity, *R* is the impeller radius, Γ<sup>0</sup> is the digester walls and Γ<sup>1</sup> 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 *τ* (Pa) as a function of the shear rate *γ*\_ (s�<sup>1</sup> ). It is calibrated by two parameters: the consistency index *k* (Pa.s<sup>n</sup> ) and the flow index *n*. It is expressed as follows:

$$
\pi = k \dot{\mathbf{y}}^n \tag{12}
$$

The viscosity is expressed as follows:

$$
\eta = k \dot{\boldsymbol{\gamma}}^{n-1} \tag{13}
$$

The components of the tensor of the shear rate are:

$$
\dot{\gamma}\_{ij} = \left(\frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i}\right) \tag{14}
$$

*CFD Simulations in Mechanically Stirred Tank and Flow Field Analysis: Application… DOI: http://dx.doi.org/10.5772/intechopen.93926*

$$\mathbf{S}\_{\vec{\eta}} = \eta \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right) = k \dot{\boldsymbol{\gamma}}^{n-1} \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right) \tag{15}$$

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]:

$$\text{Re} = \frac{\rho \text{Nd}^2}{\mu} \tag{16}$$

where *N* is the impeller speed (rev.s�<sup>1</sup> or rps), *ρ* is the fluid density (kg.m�<sup>3</sup> ), *d* is the impeller diameter (m) and *μ* is the dynamic viscosity (Pa.s).

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

$$\text{Re}\_{\mathcal{g}} = \frac{\rho U\_{\infty}^{2-n} d^n}{K \left(0.75 + \frac{0.25}{n}\right)^n 8^{n-1}} \tag{17}$$

where *U*<sup>∞</sup> is the average velocity of the fluid (m.s�<sup>1</sup> ), *k* is the consistency index (Pa.s<sup>n</sup> ), *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\pi NT\tag{18}$$

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

$$T = \left( \int\_{\mathcal{S}} \left( \vec{r} \* (\overline{\mathbf{r}} \cdot \hat{n}) \right) dS \right) \cdot \hat{a} \tag{19}$$

where *S* represents the surfaces including the rotating parts, *τ* Is the total stress tensor, *n*^ is a unit vector normal to the surface, *r* ! is the position vector and *α*^ is a unit vector parallel to the rotation axis.

The mixing energy level (MEL) (W.m�<sup>3</sup> ) is calculated by [21]:

$$\text{MEL} = \frac{P}{V} \tag{20}$$

where *P* is the power consumption (W) and *V* is the working volume (m<sup>3</sup> ).

The flow rate through the moving zone *Q* (m<sup>3</sup> .s�<sup>1</sup> ) 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 *NP* equal to 0.3 and *NQ* equal to 0.56 [11, 22]. The power and flow numbers are respectively expressed by:

$$N\_P = \frac{P}{\rho \mathbf{N}^3 d^5} \tag{21}$$

$$N\_Q = \frac{Q}{Nd^3} \tag{22}$$

The spatially average characteristic velocity gradient *G* (s�<sup>1</sup> ) 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]:

$$\mathbf{G} = \left(\frac{\overline{\boldsymbol{\varepsilon}}}{\nu}\right)^{1/2} \tag{23}$$

where *ν* is the kinematic viscosity and *ε* (m<sup>2</sup> .s�<sup>3</sup> ) is the global average turbulent energy dissipation rate calculated as follow [22, 23]:

$$
\overline{\varepsilon} = \left(\frac{N\_P N^3 d^5}{V}\right) \tag{24}
$$

The circulation time *tc* is expressed by [22, 24]:

$$t\_c = \frac{V}{N\_Q N d^3} = \frac{V}{Q} \tag{25}$$

The vorticity *ξ* (s�<sup>1</sup> ) 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 \tag{26}
$$

The helicity is expressed as:

$$H = U \cdot \xi = U \cdot (\nabla \times U) \tag{27}$$

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

$$a = \cos^{-1}\left(\frac{H}{|\xi||U|}\right) \tag{28}$$

*CFD Simulations in Mechanically Stirred Tank and Flow Field Analysis: Application… DOI: http://dx.doi.org/10.5772/intechopen.93926*
