CFD Simulation Scientific Paper (By:Jorge Rodriguez, Yong Sheng Khoo)

MSWord version

Title: CFD Analysis of a Flocculation Tank for Sustainable Drinking Water Treatment

Abstract

Flocculation is an important process used by AguaClara to treat water. The process involves particle collisions and agglomeration to form flocs. Computational Fluid Dynamics was used to better understand the fluid dynamics in the reactor. The standard K-ε model was used for every simulation model. The pressure coefficient drop over one baffle turn is 3.75, which agrees with literature estimates. After a clearance height of one baffle width or greater, the pressure coefficient drop and the maximum velocity become approximately constant. Most of the energy dissipation occurs in the region after the turn, over a distance of two baffle widths. An area of flow recirculation occurs near the center wall immediately after the turn. The pressure drop is not sensitive to the Reynolds number for a large range of inlet velocities. Better understanding of the flocculation dynamics will enable optimized particle agglomeration and break-up.

1. Introduction

Flocculation is the process by which particles collide and agglomerate. Past research has shown that shear gradients play an important role during flocculation. This process was simulated using Computational Fluid Dynamics (CFD). The main task of this research is to find the optimum strain rate in the reactor to influence particle collision. Gambit and FLUENT were utilized to model one baffle turn. Gambit was used to create the geometry of the flocculator, and to generate the mesh. FLUENT was used set up the boundary conditions and to obtain the results.

2. Methodology

The flocculation tank used by AguaClara involves 180 deg turns over few dozens baffles. To save computational effort, only one 180 deg turn was modeled. The first step was to set up the geometry of the turn. To enable future comparison with the experimental data, the geometry mimics the pilot plant. The modeling approach was to create the geometry, mesh it, set boundary conditions, and solve it using FLUENT.

2.1 Creating Geometry

Figure 1. Geometry of Flocculator

The design parameters used are:

Height: 1 m
Clearance: 0.15 m
Baffle width: 0.1 m

With these parameters, the mesh be built.

2.2 Setting up Mesh

Figure 2. Meshing Parameters (Click on the figure to see the original size)

Figure 2 shows the meshing parameters that were used. The boundary layers were first established at all the wall surfaces. They were set such that the solution would provide a result of y+ less than 5. After that, the edges were meshed as to provide higher resolution near the turn. In the final step, all the faces were then meshed.

Figure 3. Mesh of the Model (Click on figure for original size)

Figure 3 shows the mesh of the flocculator model. As can be seen, the mesh is finer near the turn and at the walls. The next step was to set up the boundary conditions of the system.

2.3 Boundary Conditions

Figure 4. Boundary Conditions

Figure 4 shows the boundary conditions used for modeling. For a flocculator, there is one in-flow and one out-flow boundary conditions. Since inlet velocity inlet was known from the experimental data, the inlet was set to the Velocity Inlet type boundary condition. The outlet was set to Pressure Outlet boundary condition type, equal to the atmospheric pressure.

The 2-Dimensional mesh was then exported to FLUENT for analysis.

2.4 Solve using FLUENT

At this stage, the Standard k-ε turbulence model was set up. Water was defined as the working material from the FLUENT database. The discretization method for the momentum, turbulent kinetic energy, and turbulent dissipation rate were set to the 'Second Order Upwind' scheme to obtain a 'Second Order Accurate' solution. The boundary conditions were set according to the values shown in the table 1. The solution was obtained by iterating until the residuals converged to 10e-6. Results were then analyzed and plotted.

Table 1. Boundary Conditions

2.5 Mesh Sensitivity Analysis

The effect of the number of mesh elements on the result was carried out. Coarse, medium and fine meshes were created and the pressure coefficient drop was compared. This analysis will provide confidence on the accuracy of certain mesh. Table 2 below shows the summary of the 3 meshes created to perform this analysis. Please refer back to figure 2 for corresponding meshing parameters.

Table 2. Mesh Meshing Parameters for Coarse, Medium and Fine Meshes

2.6 Effect of Reynolds Number

Since the inlet flow rate can vary substantially across AguaClara plants the effects of inlet Reynolds number on pressure coefficient drop was also examined.

2.7 Parameterization

At the later stage of project, the effect of geometry parameters on the results was analyzed. Different clearance heights were used for analyzing pressure drops and maximum velocities. A parameterization technique was used to automatically create a mesh given the parameters of the geometry. Using this method, the clearance height, baffle width and baffle length were easily adjusted. The Gambit journal file is included in the Appendix A.

2.8 Comparing Turbulence Model

Pressure coefficient drops were compared for the Standard K-ε, K-ε Realizable and K-ω turbulence models.

3. Results and Discussion

The results considered were plots of the velocity vectors, pressure coefficient contours, contours of strain rate and contours of turbulence dissipation rate.

Figure 5. Velocity Vectors (Click on figure for original size)

The velocity vector plot shown above depicts the water velocity throughout the flocculator. As can be seen, there is high velocity at the outer side of the turn and recirculation near the center of the wall. Furthermore, there is a region of stagnant water at the bottom of the flocculator.

Figure 6. Contours of Stream Function (Click on figure for original size)

The contours of stream function shown above tell us how particles of fluid travel in the flocculator. There is an enclosed streamline at the inner side of the turn. This means there is recirculating fluid 'trapped' in that region. The expected contraction coefficient for a 90 degree turn is 0.61. For a 180 degree turn this contraction coefficient is squared to obtain a value of 0.37. From the stream function measured at the narrowest point of the contraction the contraction coefficient is 0.44.

Figure 7. Contours of Pressure Coefficient (Click on figure for original size)

Figure 7 shows that most of the pressure coefficient drop occurs around the bend. The pressure coefficient drop is about 3.75 across the bend. This is in excellent agreement with literature estimates. The pressure coefficient drop from the expansion can be calculated analytically and is Add this equation and calculation. The difference between the measured pressure coefficient drop and the analytical solution based on the expansion is the additional pressure drop caused by shear on the walls.

Figure 8. Contours of Strain Rate (Click on figure for original size)

Contours of strain rate show high velocity gradients around the turn. There are also high strain rates at the boundary layers near the wall. It is postulated that the flocculation collision rate is directly proportional to the strain rate. The high strain rate on the wall near the contraction is of particular interest because of the extremely high values that could be responsible for floc breakup.

Figure 9. Contours of Turbulent Dissipation Rate (Click on figure for original size)

Contours of turbulent dissipation rate show a similar trend as the contours of strain rate right after the turn. The region of highest turbulence dissipation occurs after the turn. As can be seen from figure 9, this region is roughly twice the length of baffle spacing.

An analytical estimate of the energy dissipation rate in the dissipation cell can be obtained. Use the following equations:



Using the equation above to estimate the energy dissipation rate in the expansion zone we obtain an average value of ??. This is based on the assumption that the length of the dissipation region is approximately twice (Π _cell = 2) the distance between the baffles.

Calculate the value of Gθ for each cell in the computational grid and sum this up over the domain containing the dissipation zone to get a more accurate measure of the actual Gθ for a 180 degree bend. Compare that with the estimate that the AguaClara team is currently using.

Figure 10. Wall Yplus

Figure 10 shows that the y+ values. According to FLUENT documentation "the mesh should be made either coarse or fine enough to prevent the wall-adjacent cells from being placed in the buffer layer (y+ = 5~30)". Since the y+ from the model was consistently less than five (in the viscous sublayer) the turbulence near the walls was resolved properly.

Figure 11. Mesh Sensitivity Analysis

Figure 11 shows the pressure coefficient drop over one turn for different mesh densities. In general, finer meshes provide more accurate results. However, as the mesh was refined the pressure drop remained constant as can be seen in figure 11. Hence, it was concluded that results were not sensitive to mesh density and the coarse mesh was sufficient.

Figure 12. Reynolds Number Effect on Pressure Coefficient Drop

The effect of the Reynolds number on the pressure coefficient drop was analyzed. This was done by changing the inlet velocity which initially produces a Reynolds number of 10,000. From figure 12, it can be seen that the value of the pressure coefficient drop has a small change when compared to big changes in Reynolds number. In other words, the pressure coefficient drop is not sensitive to the Reynolds number at the inlet. This implies that the design of the flocculator should not be altered by the inlet flow rate. This is to say that one flocculator design can be used for different flow rates. You are right that the pressure coefficient is insensitive to Reynolds number. But that does not mean that a single flocculator can handle a wide range of flows. The energy dissipation rate is still proportional to the velocity cubed.

Figure 13. Clearance Height Effect on Pressure Coefficient Drop and Maximum Velocity

The effect of the clearance height on the pressure coefficient drop was also analyzed. It can be seen that the pressure coefficient drop is independent of the change in clearance height after the clearance height is greater than a critical value. Figure 13 shows that after a critical value of 1, the pressure coefficient drop is constant. This phenomenon can be explained by looking at figure 14 below. Figure 14 shows the turbulent dissipation rate for clearance heights of 1b and 1.5b. It can be observed that the maximum dissipation rate is equal for both reactors. A recirculation zone begins to form at the bottom of the reactor when the clearance height is larger than 1b. A similar argument can be made for the values of maximum velocity. A clearance height less than 1b results in a higher pressure coefficient drop as it creates an additional constriction in the flow increasing expansion losses. It is therefore recommended for the design team that the clearance height be at least the same as the baffle spacing. The correlation between pressure coefficient drop and maximum velocity should also be noted. Add the theoretical connection using the pressure drop in an expansion. The conclusion that there is no need to make the clearance height larger than 1b is noteworthy. This will make it possible to handle higher flow rates with the same channel width since we are dealing with a constraint that requires neighboring baffles to overlap in the center of the reactor to prevent short circuiting.

Figure 14. Comparison of Turbulent Dissipation Rate for Clearance height of 1b and 1.5b

Figure 14 above further validates that results are not sensitive to the change in clearance height. Contours of turbulence dissipation rate for clearance heights of 1b and 1.5b were compared. These results show that the region of active turbulent dissipation is the same for both reactors, about two times the length of baffle spacing.

Figure 14. Effect of Turbulence Model on Pressure Coefficient Drop

Pressure coefficient drop was least in the K-ε model and most in the K-ω model.

Figure 15. Contours of Velocity Magnitude for Different Turbulence Models

As figure 15 depicts the standard K-ε has the smallest region of high and low velocities, shown in red and blue respectively. This explains why it has the lowest pressure coefficient drop. The K-ω model has the biggest region of high and low velocities therefore having larger pressure coefficient drops.

However, as it can be observed from the transparent demo plant, the recirculation area is only the length of one or two baffle widths. This is contrary to the excessively large blue/red regions predicted by the k-ε realizable, and k-w models. Therefore it was concluded that the standard k-ε model best simulates the turn.

Figure 16. Contours of Stream Function for Different Turbulence Models

Figure 16 clearly shows the recirculation region for the three different turbulence models. Since standard K-ε model best represent the flow features seen in the demo plant, it is concluded that standard K-e model best simulate the turn.

I am uncomfortable with the method of selecting the turbulence model. The demonstration plant is in the laminar flow region and thus it may not be a good reference point for full scale turbulent reactors. The difference between the various models is quite significant especially when you consider that we are most interested in Gθ. The models that predict a large wake behind the expansion will also predict lower average values of the energy dissipation rate and that will correlate with a larger Gθ. Thus we need a better basis for choosing the model that we will use for further research. I suggest further research to determine which models do a better job of modeling the evolution of a turbulent jet. The region where the models diverge significantly is in their predictions about the length of the zone influenced by the jet.

4. Conclusions

1. An area of recirculation occurs near the center wall immediately after the turn
2. Pressure coefficient drop over one baffle turn is 3.75
3. After a clearance height of one baffle spacing or greater, the pressure coefficient drop and the maximum velocity becomes constant
4. Most of the energy dissipation occurs in the region after the turn over a distance of about two baffle spacings
5. The pressure coefficient drop is insensitive to the Reynolds number for a large range of inlet velocities
6. A mesh with 20,000 mesh elements is sufficient to obtain accurate results
7. Different turbulence model resulted in fairly different results
8. The standard k-e model best simulates the turn

6. Future Research

1. Measure the Gθ value from FLUENT
2. Analyze the region of very high energy dissipation along the wall near the contraction and explore methods to reduce the energy dissipation in this zone
3. Devise improved methods of creating more uniform energy dissipation in the reactor to enhance flocculation efficiency
4. Analyze the effect of changing the ratio of b/h all the way to the extreme where the baffles don't overlap in the center of the reactor. This is important to learn what will happen as we design for larger flow rates. Larger flows require the spacing between baffles to increase. If we hold the reactor depth constant that will increase the ratio of b/h.
5. Better understanding of different turbulence models
6. Model droplet collision - breakup

7. Acknowledgements

During the course of this project, we received invaluable advice and direction from Prof. Monroe Weber-Shirk, whose leadership of the AguaClara team project spurred this research and technical investigation. All questions we had were directed towards Dr. Rajesh Bhaskaran, his expertise and academic guidance kept us on schedule and on task. We would also like to thank Prof. Brian Kirby for elucidating some technical aspects we encountered. For their help, we are grateful.

References

On the collision of drops in turbulent clouds P. G. Saffman and J. S. Turner, Journal of Fluid Mechanics Digital Archive, Volume 1, Issue 01, May 1956, pp 16-30 doi: 10.1017/S0022112056000020, Published online by Cambridge University Press 28 Mar 2006

Appendix A

/ Journal File for GAMBIT 2.4.6, Database 2.4.4, ntx86 SP2007051421
/ Identifier "Clearance 1.5W"
/ File opened for write Fri Apr 18 10:07:42 2008.
undo begingroup

$clearance = 0.1

coordinate modify "c_sys.1" xyplane xaxis add 0 AND 0.1 AND 0.2 reset snap \
lines
coordinate modify "c_sys.1" xyplane yaxis reset snap lines
window modify coordinates "c_sys.1" xyplane grid
undo endgroup
undo begingroup
coordinate modify "c_sys.1" xyplane xaxis add 0 AND 0.1 AND 0.2 reset snap \
lines
coordinate modify "c_sys.1" xyplane yaxis add 0 AND 0.05 AND 0.1 AND 0.15 AND \
0.2 AND 0.25 AND 0.3 AND 0.35 AND 0.4 AND 0.45 AND 0.5 AND 0.55 AND 0.6 AND \
0.65 AND 0.7 AND 0.75 AND 0.8 AND 0.85 AND 0.9 AND 0.95 AND 1 reset snap \
lines
window modify coordinates "c_sys.1" xyplane grid
undo endgroup
vertex create coordinates 0 0 0
vertex create coordinates 0.1 0 0
vertex create coordinates 0.2 0 0
vertex create coordinates 0.2 $clearance 0
vertex create coordinates 0.1 $clearance 0
vertex create coordinates 0 $clearance 0
vertex create coordinates 0.2 1 0
vertex create coordinates 0.1 1 0
vertex create coordinates 0 1 0
edge create straight "vertex.4" "vertex.7"
edge create straight "vertex.5" "vertex.8"
edge create straight "vertex.6" "vertex.9"
edge create straight "vertex.8" "vertex.9"
edge create straight "vertex.7" "vertex.8"
edge create straight "vertex.3" "vertex.4"
edge create straight "vertex.2" "vertex.5"
edge create straight "vertex.1" "vertex.6"
edge create straight "vertex.5" "vertex.6"
edge create straight "vertex.4" "vertex.5"
edge create straight "vertex.2" "vertex.3"
edge create straight "vertex.1" "vertex.2"
undo begingroup
coordinate modify "c_sys.1" xyplane xaxis add 0 AND 0.1 AND 0.2 reset snap \
lines
coordinate modify "c_sys.1" xyplane yaxis add 0 AND 0.05 AND 0.1 AND 0.15 AND \
0.2 AND 0.25 AND 0.3 AND 0.35 AND 0.4 AND 0.45 AND 0.5 AND 0.55 AND 0.6 AND \
0.65 AND 0.7 AND 0.75 AND 0.8 AND 0.85 AND 0.9 AND 0.95 AND 1 reset snap \
lines
window modify coordinates "c_sys.1" xyplane nogrid
undo endgroup
face create wireframe "edge.4" "edge.2" "edge.9" "edge.3" real
face create wireframe "edge.5" "edge.1" "edge.10" "edge.2" real
face create wireframe "edge.6" "edge.11" "edge.7" "edge.10" real
face create wireframe "edge.12" "edge.7" "edge.9" "edge.8" real
undo begingroup
/ERROR occurred in the next command!
blayer create first 0 growth 1.2 total 0 rows 4 transition 1 trows 0 uniform
undo endgroup
undo begingroup
blayer create first 0.0003 growth 1.25 total 0.0077407 rows 9 transition 1 \
trows 0 uniform
blayer attach "b_layer.1" face "face.2" "face.3" "face.3" "face.4" "face.1" \
"face.2" "face.1" "face.4" "face.4" "face.3" edge "edge.1" "edge.6" \
"edge.7" "edge.7" "edge.2" "edge.2" "edge.3" "edge.8" "edge.12" "edge.11" \
add
undo endgroup
undo begingroup
edge picklink "edge.3" "edge.2" "edge.1"
edge mesh "edge.1" "edge.2" "edge.3" successive ratio1 1.01 size 0.005
undo endgroup
undo begingroup
edge picklink "edge.6" "edge.7" "edge.8"
edge mesh "edge.8" "edge.7" "edge.6" successive ratio1 1 size 0.002
undo endgroup
undo begingroup
edge picklink "edge.4" "edge.5" "edge.10" "edge.9" "edge.11" "edge.12"
edge mesh "edge.12" "edge.11" "edge.9" "edge.10" "edge.5" "edge.4" successive \
ratio1 1 size 0.002
undo endgroup
face mesh "face.2" map size 1
face mesh "face.1" map size 1
face mesh "face.3" map size 1
face mesh "face.4" map size 1
physics create "Inlet" btype "VELOCITY_INLET" edge "edge.4"
physics create "Outlet" btype "PRESSURE_OUTLET" edge "edge.5"
physics create "Left_Wall" btype "WALL" edge "edge.3" "edge.8"
physics create "Middle_Wall" btype "WALL" edge "edge.2"
physics create "Right_Wall" btype "WALL" edge "edge.1" "edge.6"
physics create "Bottom_Wall" btype "WALL" edge "edge.12" "edge.11"
export fluent5 "Clearance 1.5W.msh" nozval
save