Abstract

This work shows the achievements obtained in the numerical optimization of mixing operations in water treatment plants thanks to the development of a Finite Element hydrodynamic code that evaluates and optimises the hydrodynamics of such processes. Some water treatment units present important deficiencies in their design. The numerical optimization of the flow manages to avoid the high expenses of the trial and error processes involved in building, installing and removing the mixing mechanisms, channels and basins, and those derived from the need of halting the water treatment processes. As a result, an optimum design of the treatment plant is obtained at a low cost.

Keywords

Agitation, Finite Element Method, Mixing, Navier-Stokes, Optimization, Water Treatment, Multidimensional Constrained Direct Optimization

Share and Cite:

1. Introduction

Some of the most challenging problems of the 21st century are those related to environmental issues. Among those, the efficiency in the energy consumption and an adequate management of the water resources are undoubtedly key aspects. The aim of the present research attempts to address the search for an energy-efficient optimization of a water treatment process through the programming of a Finite Element code that not only evaluates the hydrodynamic behaviour of the water in the water treatment basin, but also incorporates a mathematical optimization of the involved parameters in order to obtain their most efficient performance.

Some previous studies have addressed this topic by several researchers. In particular, a previous study by the authors was published in [1] regarding the advances that could be introduced into mixing processes with the help of a hydrodynamic evaluation of the flow. Nonetheless, this previous work would only evaluate the improvements to be obtained on a trial-and-error basis, through a sequential resolution of the FE code developed in that work at different scenarios. The present work shows the development of a MATLAB code that not only solves the water flow in an efficient way through the Finite Element Method, but also achieves a mathematical optimization of a particular mixing process in order to obtain the optimal values for the geometrical and physical parameters defined as the key values defined by the user.

The development of the code involves not only the well-known difficulties associated with the resolution of the viscous incompressible flow, but also those that show up when incorporating optimization tools. The code used in the numerical evaluation of the hydrodynamics of the treatment tanks has been a new release of the HYDRAFEM code (stable HYDRodynAmic Finite Element Model), first developed in [2]. Still some new features have been introduced together with a new MATLAB implementation, as it will be further explained in Chapter 2.

Once the hydrodynamics of the flow have been obtained for the present conditions, the main design parameters are defined and an optimization tool is programmed in order to obtain the best choice for those parameters in order to achieve an optimum behaviour of the mixing tank. The optimization procedure has been based upon a multidimensional constrained direct optimization method with a parabolic interpolation for each of the optimization variables as it will be described in Chapter 3.

In this work, the hydrodynamic behaviour of certain neutralizing tanks with an inflow of some 5 m³/s will be optimized. Each of the tanks is equipped with three mixing impellers. The plan dimensions of the tanks are 20.5 × 7 m² with a water depth of some 4.2 m. These dimensions match up those of the neutralizing units used in the As Pontes water treatment plant. This treatment plant owned by the firm ENDESA, was devoted to the treatment of the water outflow obtained from the opencast coal mine of As Pontes, feeding the homonymous power station. Neutralizing units of similar characteristics can be found in many water treatment plants all over the world. The particulars about the mixing units can be found in Chapter 4. The aim of this work will be to optimize the efficiency of the neutralizing processes. The results regarding this particular case can be seen in Chapter 5.

2. Finite Element Resolution of the Fluid Flow

A new code for the optimization of the hydrodynamic behaviour of the flow in water basins has been developed. The code is based on a new version of the hydrodynamic model, HYDRAFEM, developed by the first author [3] in which apart from a newly developed optimizing tool, some other new features have been included.

HYDRAFEM makes an approximation of the behaviour of flows with a moderate Reynolds number, such as those that take place at water treatment plants. There were some previous works that evaluate the flow in water treatment units by using the Stokes simplification (“creeping flow”) or the potential flow equations [4]. This type of simplifications allows for the obtaining of an approximate solution of the flow. Some other interesting approaches for the optimization of water treatment plants performances are considered at [5] and [6]. In order to obtain the real behaviour of the flow, the solution of the full Navier-Stokes equations is needed, and the numerical resolution of these equations involves tough stability problems which do not show up for the former simplified formulae. The code has previously been applied to several problems related to Environmental Engineering with very interesting results [7]. A work was done in 2007 to improve the behaviour of mixing operations in water treatment plants [1] on a trial and error basis. A newly developed numerical optimization tool has been incorporated to the code and its achievements have been checked on the same mixing operations used in [1] as a practical case. Details regarding the optimization will be further explained in section 3 to later present the particulars about the hydrodynamic model. A new program has been coded in Matlab for this application that differs from the previous HYDRAFEM code, previously programmed in FORTRAN. The particulars are as follows.

The equations that govern the viscous incompressible flow are the so-called Navier-Stokes equations. Made up by the dynamic and the continuity equations, they are named after their discoverer, the French civil engineer Claude-Louis Navier (1785-1836), who in 1821 formulated the equations that rule the incompressible flow. Using the indicial notation we can express the Navier-Stokes equations as:

$u_{i,t}+u_j{u}_{i,j}=-\frac{1}{\rho }{p}_{,i}+\nu u_{i,jj}+f_i$

$u_{i,i}=0$ (1)

in the domain of integration $\Omega$ , together with the boundary and initial conditions:

${u_i]}_{{\Gamma }_1}=b_i,{{\sigma }_{ij}{n}_j]}_{{\Gamma }_2}=t_i$

$u_i\left(x_j,0\right)=u_{i0}\left(x_j\right)$ (2)

where u_i is the velocity, p is the pressure, f_i is the body force per unit mass, $\rho$ is density, $\nu$ is the cinematic viscosity, ${\Gamma }_1$ and ${\Gamma }_2$ are two non-overlapping subsets of the piecewise smooth domain boundary $\Gamma$ , b_i is the velocity vector prescribed in ${\Gamma }_1$ , t_i is the traction vector prescribed on ${\Gamma }_2$ , ${\sigma }_{ij}$ is the stress along the boundary ${\Gamma }_2$ , and n_j is the outward unit vector normal to ${\Gamma }_2$ .

The unknowns in these equations are both the velocity and the pressure. The most straight forward way of dealing with these two kinds of variables when using the Finite Element Method is to use the so-called mixed methods, in which both of them are considered as the unknowns of the resulting system of equations. The corresponding weighted formulation is obtained from the application of the weighted residuals method on both the dynamic and continuity equations. A Galerkin formulation is used in which the test functions are taken equal to the trial functions, in terms of which the unknowns are interpolated. A zero-equation turbulent model is used by increasing the magnitude of the kinematic viscosity to account for the value of the turbulent viscosity.

$\begin{array}{c}{\displaystyle \underset{{\Omega }_h}{\int }{w}_i^h\left(u_{i,t}^h+u_j^h{u}_{i,t}^h-f_i^h\right)}d\Omega +v{\displaystyle \underset{{\Omega }_h}{\int }{w}_{i,j}^h{u}_{i,j}^h}d\Omega -\frac{1}{\rho }{\displaystyle \underset{{\Omega }_h}{\int }{w}_{i,i}^h}pd\Omega -{\displaystyle {\int }_{{\Gamma }_2}{t}_i^h{w}_i^{h}d{\Gamma }_2}\\ +{\displaystyle \sum _e{\displaystyle \underset{{\Omega }_e}{\int }{p}_i^h\left(u_{i,t}^h+u_j^h{u}_{i,t}^h-v{u}_{i,jj}^h+\frac{1}{\rho }{p}_i^h-f_i^h\right)d\Omega }}=0\end{array}$ ;

${\displaystyle \underset{{\Omega }_e}{\int }{q}^h{u}_{i,i}^h}d\Omega =0$ (3)

where w_i are the functions for the dynamic equation and q are the weighting functions for the continuity equation. p_i will be taken equal to w_i in the Galerkin formulation. The h superscript stands for the discretization being carried out in our formulation, which will be made in terms of a P2P1 basic element (see Figure 1).

Figure 1. P2P1 basic element.

This formulation should verify in addition to the governing equations, some consistency conditions, the most restrictive of which is the LBB (Ladyzhenskaya-Babuska-Brezzi), or divergence-stability condition. An inappropriate choice of this basic element may fail to verify this condition and prevent the whole algorithm from converging. The P2P1 basic element, that is quadratic velocity, linear pressure, has been shown to produce stable and efficient solutions [8]-[11]. The Q1Q0 basic elements used in the previous formulation have been avoided in order to prevent problems such as the pressure checkerboard.

Once the elementary matrices are evaluated and assembled, the integral Equation (3) can be expressed in matrix notation as:

$M_{\upsilon }\frac{\partial \underset{\_}{\upsilon }}{\partial t}+C_{\upsilon }\left(u,v\right)\underset{\_}{\upsilon }+\nu A_{\upsilon }\underset{\_}{\upsilon }-Bp=f$

$B^T\underset{\_}{\upsilon }=0$ (4)

where $M_{\upsilon }$ is the mass matrix, $C_{\upsilon }\left(u,v\right)$ is the convective matrix, $A_{\upsilon }$ is the viscous matrix, B is the pressure matrix, f is the external forces vector, p is the pressure vector, u is the velocity vector in the x direction, v is the velocity vector in the y direction and $\underset{\_}{\upsilon }$ is the velocity vector.

For the time derivatives included in the first term of the first equation in (4), a Finite Difference backward differencing approach will be used in order to transform our partial differential system into an algebraic one. The convective term that appears in (4) is a non-linear velocity-dependent function. This term should be eliminated in order to transform system (4) into a linear system of equations. The numerical scheme to be used for this transformation is a so-called successive approximation method that can be mathematically expressed as:

${\displaystyle \underset{{\Omega }_h}{\int }{\overline{w}}_i^{}{u}_j^{}{u}_{i,j}^{}}d\Omega \approx {\displaystyle \underset{{\Omega }_h}{\int }{\overline{w}}_i^{}{u}_j^{n-1}{u}_{i,j}^n}d\Omega$ (5)

3. Numerical Optimization

Once the hydrodynamics of the flow have been solved, the optimization procedures are introduced. The aim would be to obtain an n dimensional design vector, y, that minimizes or maximizes a certain objective function f(y). The variables in y will be subject to certain constraints that can be written as

$d_i\left(y\right)\le ai,i=1,2,\cdots ,m$

$e_i\left(y\right)\le bi,i=1,2,\cdots ,p$ (6)

The optimization is performed using a direct (derivative not being assessed) multidimensional univariate search method, that is, each of the variables will be changed one at a time to improve the approximation while the other variables are held constant, thus reducing to a sequential one-dimensional search. For each approximation, a parabolic interpolation of each of the variables will be used (see Figure 2).

Figure 2. Univariate multidimensional optimization with parabolic interpolation.

The parabolic interpolation takes advantage of the fact that a second-order polynomial usually provides a good approximation to the shape of the objective function near an optimum. The constraints to be imposed will normally allow for the bracketing of the design variable within a certain interval. A third point in-between (within a golden ratio) will be selected and the maximum of the three-point parabola will be determined as

$x_3=\frac{f\left(x_0\right)\left(x_1^2-x_2^2\right)+f\left(x_1\right)\left(x_2^2-x_0^2\right)+f\left(x_2\right)\left(x_0^2-x_1^2\right)}{2f\left(x_0\right)\left(x_1-x_2\right)+2f\left(x_1\right)\left(x_2-x_0\right)+2f\left(x_2\right)\left(x_0-x_1\right)}$ (7)

Depending on the relative position of the parabola maximum within the interval, one of the limits will be dropped and the process will be repeated until convergence.

The objective function to be considered for this problem will be the minimum typical deviation of the velocity module x in the basin as it will be further explained in the chapter to follow. The standard deviation can be written as the square root of the variance of n values with $\overline{x}$ mean as

$f\left(y\right)=\sigma \left(y\right)=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(x_i\left(y\right)-\overline{x}\left(y\right)\right)}^2}}{n}}$ (8)

4. Study Case: Agitation Processes

The hydrodynamic behaviour of the water treatment units installed in many plants shows very important deficiencies. The adequate operation of these treatment units is closely related to the proper hydrodynamic behaviour of the flow. The physicochemical processes taking place in these units should be carried out within the frame of correct global hydrodynamics, which should be adapted to the characteristics that feature the treatment plant. The numerical optimization of the hydrodynamic behaviour of the basin will allow for an efficient energy consumption of the water treatment plant leading to an important decrease in operating costs.

The study case considered for this work is a mixing basin in a wastewater treatment plant, where the mixing tank is mainly used for neutralization purposes. The addition of lime reduces corrosiveness from acid water raising the pH. The work is based on a previous study with no optimization involved and the same dimension and parameters will be assumed for a start.

The neutralizing unit considered in the present work will be, as in [1], consisting of two sets of three tanks each. The inlet discharge at each of the inlets is estimated in some 5 m³/s, resulting in a mean velocity at the inlet of about 0.75 m/s. The water volume is of about 600 m³, and the depth has been taken as 4.5 m. The axial symmetries in the tank allow us to use a 2D evaluation of the vertical section of the flow as very appropriate for our purposes.

In order to achieve an adequate neutralization of the effluent, contact among particles should take place at every single part of the tank. The physicochemical mechanisms that permit the mixing processes are convection, turbulent diffusion and molecular diffusion [12]. The agitation processes should also prevent the settling of the suspended solids that are usually present in the water effluent. Therefore, the velocity should be kept to a minimum while an efficient agitation in terms of energy consumption is maintained. This energy-efficient agitation will keep the solids in suspension and will remove solids from the bottom while letting the migration and the mixing processes of the liquid phase take place within the effluent [13].

An axial pattern mixing element is the most adequate to achieve this energetic and efficient agitation. The radial (b) and tangential (c) flow patterns (see Figure 3) are those that follow respectively the direction of the radius and the tangent of the circle plotted by the propeller in its movement. In the axial flow pattern the main flow direction is that of the agitator axis and can be propelling or sucking depending on the direction of the spin [14]. An axial mixing element with a diameter D, produces a stream of water whose radius r_c is bigger than D/2 (see Figure 4). Due to the position of the bottom of the tank, the streamlines bend slightly outwards.

Figure 3. Flow patterns. (a) Axial, (b) Radial, (c) Tangential.

Figure 4. Axial velocity field and streamlines.

As a consequence, a recirculation region is generated at both sides of the agitator and a “dead cone” shows up underneath the agitator [15]. The size and magnitude of the recirculation streamlines will depend on the inflow discharge, the volume of water provided by the agitator, the position of the propellers and the geometry of the tank.

Propellers with slender blades, such as the SABRE® (Dosapro Milton Roy) blades, as the one shown in Figure 3, are well suited to the homogenization of an effluent with suspended solids. They produce an efficient flow with scant energy consumption. The so-produced flow is mainly axial, having its radial and tangent components reduced to a minimum. The resulting turbulent percentage is from very low to moderate, contributing in this way to a good energy efficiency. This is why this kind of mixing element is very adequate for applications in which high velocities at a low cost are required, such as neutralizing processes of big volumes with suspended solids. With the tangential and radial velocities being kept to a minimum, the effect of the propellers can be characterized in terms of a given constant vertical velocity.

The limiting parameter governing the agitation is the pumping volume of each of the agitators. The neutralizing process is solved by laying out three agitators in each of the tanks. The calculated pumping volume has been evaluated in some 19,000 m³/h, which, when used with propellers with a diameter of 2.5 m results in an axial flow velocity of about 1 m/s. So as to reach these working parameters, 11 kW motors with 1500 rpm and a differential of 60 rpm have to be used, with an elevation from the bottom of around 1.2 m.

5. Numerical Results

The numerical evaluation was performed using a P2P1 triangular mesh with 441 nodes and 192 elements as shown in Figure 5. For the present working conditions an inlet discharge of 5 m³/s has been considered on the left-hand side. The results obtained for the velocity field have been those shown in the figures to follow.

In Figure 6, the position of four velocity profiles placed along the depth of the tank can also be observed. These velocity profiles provide information about the velocity of the flow at depths a = 0.6, b = 1.8, c = 2.4 and d = 3.6 m, allowing evaluation of the proper operation of the mixing processes at significant points.

Figure 5. Finite element mesh.

Figure 6. Velocity profiles.

The velocities at these profiles will also be used to assess an objective function for the optimization as the standard deviation, of the velocity at the vertical section of the tank. With the objective function being

$f\left(x\right)=\sigma \left(v_0,d\right)=\sqrt{\frac{{\displaystyle \sum _{i=a}^d{\left(x_i\left(v_0,d\right)-\overline{x}\left(v_0,d\right)\right)}^2}}{4}}$ (9)

where $x_a,x_b,x_c,x_d$ are the mean velocities along depths a, b, c, d and $\overline{x}=\frac{x_a+x_b+x_c+x_d}{4}$ .

The limiting mixing mechanism is convection, and the characteristic value that features the homogenization rate is the pumping volume provided by the agitator. This energetic agitation may cause excessive forces against the mixing elements, which could lead to blade failure. This misbehaviour will be avoided by keeping the downward velocity at the blades underneath v₀ = 1.5 m/s as an optimization constraint. A minimum value of v₀ = 0.5 m/s has also been assessed as an optimization constraint in order to ensure a minimum mixing velocity within the basin. The downward pumping discharge range has been set between 9500 and 28,500 m³/h.

As stated before, the distance of the propellers from the bottom, d, is also a key parameter and will be considered as the second optimization parameter for the bi-dimensional design vector. Although the propellers can be placed at a distance of up to one half of the total depth from the bottom, it is still recommended that the distance of the agitators from the bottom of the tank should not be greater than 1/3 of the total depth [14]. A minimum distance from the bottom should also be ensured. This way, an allowance interval in-between d = 0.9 and d = 1.5 will be considered for the variation in this second design parameter. The design parameter constraints are therefore

$\left(x\right)=\sigma \left(v_0,d\right);0.5\le v_0;0.7\le d\le 1.7$ (10)

v₀ = The optimization problem was solved and some results are obtained within less than one minute of CPU time on a personal computer. The optimized values for the conditions considered have been a distance d of 0.9 m from the bottom and a pumping downward discharge of 28,500 m³/h (1.5 m/s). The optimized mesh with the new position for the mixer conditions then became the one shown in Figure 7.

Figure 7. Optimized mesh.

The results for velocities, streamlines and velocity profiles can be seen in Figure 8.

Figure 8. Results for optimized velocities, streamlines and velocity profiles.

The 2D approach in the vertical plane that contains the agitators axes is the most appropriate approach as it considers the sections in which the main and most critical flows are taking place with regard to the flow-behavior problems being considered. Obviously, this simplification neglects some important 3D effects, the most important of which are the shaft forces from the lateral flows being ignored. Anyway the main flow direction in the plant is the one going from the inlet to the outlet parallel to the sides of the tank, and the main flow direction in depth is the vertical one produced by the impellers in an axial flow pattern. The radial and tangential effects of the impellers are negligible when using this kind of axial agitators (see section 4), which is why considering of the vertical central section is so meaningful. Moreover, the 2D model is computationally very efficient. The flow pattern in plan elevation is expected to be mostly symmetric due to the scant influence of the tangential flow produced by the impellers, and this fact makes less important the shaft effects of the non-considered lateral “corridors” in-between the impellers and the longitudinal wall.

As can be observed, the agitation obtained is as expected. The lateral splitting of the flow takes place while a stagnant cone appears on the tank bottom, roughly concentric with the impeller axis. As can be seen in the figures the agitation pattern causes the radial splitting of the flow underneath the mixing elements. As a consequence, three primary recirculation vortices show up in which the flow velocity is relatively high. The separation of the flow governed by the agitators is asymmetric, being affected by the main direction of the flow. This increase in the velocity of the flow in a direction different from the main one is a wanted effect in the mixing processes as it allows for a further contact among the particles of the fluid, contributing in this way to their homogenisation.

Apart from the primary vortices referred above, a group of smaller secondary vortices appears in the upper part of the section, such as the one showing up at the lower right hand corner. Unlike the primary ones, these vortices are featured by very small circulation velocities and indicate the presence of stagnant regions. The existence of these recirculation regions does not contribute to the correct development of the homogenisation processes but under the present conditions they are small enough so as not be causing any major problems.

The velocity at the bottom regions (elevation velocity profiles a and b) is much more influenced by the effects of the propellers and can be easily seen at the figures. Thus, the standard deviation is decreased as the propellers are immersed into the tank, which does not necessarily increases the recirculation and the mixing time but still approaches the mean velocity values of the lower-level velocity profiles to those at the upper levels. Still minimization of the standard deviation ensures that the overall velocity across the vertical section containing the axis is as balanced as possible.

6. Conclusions

A newly developed Finite Element code has been programmed that numerically optimizes the design parameters that feature the hydrodynamic flow. The code solves the Navier-Stokes equations on a 2D basis making use of a constant turbulent viscosity that provides with the velocity and pressure field for the flow of water up to a moderate Reynolds number.

An optimization tool has been incorporated into the program in order to be able to obtain the relative optimum of a given objective function, subject to certain constraints, for a certain value of the design parameters. The optimization process has been based upon a multidimensional univariate optimization with parabolic interpolation of the design parameters.

The program has been applied to the resolution of the flow in mixing tanks such as the ones being frequently used in the water treatment industry. The main design parameters have been defined as the velocity and the position of the mixing impellers. The standard deviation of the velocity within the basin has been chosen as the objective function for driving the overall mixing process to a working optimum. The optimum values for the design parameters and the hydrodynamic variables have been finally obtained and presented.

Editorial Note

This pass focuses on English grammar, usage, clarity, and academic style. Scientific validity, equations, numerical values, and references were not independently verified.

Author Contributions

Conceptualization, P. Vellando; methodology, P. Vellando, R. Juncosa, F. Padilla and H. García-Rábade; software, P. Vellando; validation, P. Vellando; writing P. Vellando. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References