Computational Fluid Dynamics

• See the Fall 2009 Flocculator Design Case here

Flocculation Tank Simulation

Introduction

Computational Fluid Dynamics (CFD) is a tool used by AguaClara to obtain a description of the flow through a portion of the plant where particle formation and growth (flocculation) occurs. This part of the plant is referred to as the flocculator. In this section, dirty (turbid) water flows through a series of baffles that enhances turbulent mixing. Essentially, for particles (flocs) to grow and eventually settle out in the sedimentation tank, they first must collide. By increasing the level of turbulence, the flocculator is increasing the collision potential of flocs. However, if there is too much turbulence, the flocs will break up and not settle out. By running CFD simulations of the flocculator, AguaClara can analyze the parameters important to flocculation and use the resulting data when making design decisions.

Goals

The primary goal of the CFD team is to obtain an accurate solution for the flow field through the entire flocculator. To understand why the team is interested in the whole flocculator as opposed to just a section (e.g. the entrance region), it is informative to look at the constraints on the flocculation process. As mentioned in the opening paragraph, the level of turbulence must be kept between some minimum and maximum amounts to ensure floc formation and growth while preventing flocs from breaking up. Mathematically, collision potential (Ψ) is hypothesized to be proportional the turbulent energy dissipation rate (ε) to the 1/3 power. Thus, by increasing the turbulent mixing, the flocculator is increasing turbulent energy dissipation and thus collision potential as well. Conversely, high energy dissipation rates mean high strain rates. If these become too large, the particles will break apart, and effluent water will remain turbid. Further complicating flocculation is the realization from currently operating plants that flocs are settling out prior to reaching the sedimentation tank. In order to properly handle this issue, higher energy dissipation rates may be needed near the end of the flocculator to ensure that the largest flocs will break up while also creating large turbulent eddies that can pull prematurely settled flocs off the bottom of flocculator. All of these design requirements motivate the need for a flow field solution through the entirety of the flocculator.

Within the context of the primary goal comes the need to deal with all the issues that hinder finding an accurate solution, and the resolution of these issues will be considered the CFD team's secondary goal. CFD codes inherently contribute several sources of error to any numerical solution based on the structure of the computational solver. The discrete and iterative nature of CFD creates two primary sources of error for any simulation. Further, the computational model used to handle the randomness associated with turbulence and the modification of that model to deal with flow near the wall both introduce error based on their ability to properly model the problem-specific flow regime. The mitigation of the effects of these errors will be discussed in terms of the verification and validation of the flocculator model. For the purpose of CFD, verification will involve dealing with error from discretization and linearization (iteration). Validation will then look at the effects of the turbulence model compared to experimental data. In other words, verification will involve making sure the team solves the "equations right," and validation will deal with solving the "right equations." The efforts undertaken to reduce these errors will hopefully establish some standards for future AguaClara simulations.

Building Meshes in Workbench

The processes of creating geometry and meshing are fairly straightforward, so there is no need to fill this report with details on how to perform these tasks. There is, however, one trick that is worth mentioning for future years since it is non-intuitive. In order to create a uniform mesh where cell sizes are the same through the core of the flow, the thickness of each baffle needs to be infinitely thin. Unfortunately, ANSYS has only a convoluted way of doing this. Simply creating a series of lines in a sketch and converting them to line bodies will not work. For some reason, the meshing tool does not recognize line bodies. The following are the steps needed to create infinitely thin baffles.

• Step 1: Draw the rectangular region that will serve as the fluid flow domain and create the corresponding surface body
• Step 2: Create a new sketch, draw the lines that will serve as the baffles, and create line bodies from this sketch
• Step 3: Under the Tools menu at the top of the DesignModeler, select Projection
  
• Step 4: At the bottom left of the screen, a details menu will appear; for Edges, select the line bodies that will serve as baffles and click Apply
• Step 5: For Target, select the surface body that will serve as the flow domain and click Apply
  
• Step 6: Click Generate (the lightening bolt)
  

Verification and Validation

Grid Convergence Verification (Discretization Error)

Before any discussion of grid convergence, figure 1 will help clarify what is meant by baffle space # and baffle #:

Figure 1: Definition of Terms

Having defined these two terms, the discussion now turns to grid (aka mesh) size, a common source of error for computational problems. The finite volume method used in FLUENT requires that the flow domain be divided up into areas (2D) or volumes (3D) and the resulting equations be solved for each cell. Since each cell is a finite portion of the domain, there is a truncation error associated with solving the discrete equations. Decreasing cell size will help reduce this error because the discretization will asymptotically approach a point which is the same as the continuum formulation. The tradeoff comes from the computational cost associated with smaller and therefore more cells. The approach to minimizing grid error thus involves systematically decreasing the grid size by dividing the cell length and width in half and looking at the results for an otherwise constant setup.

There are essentially two methods for determining grid convergence once simulation data has been obtained. The first test is a simple qualitative look at contour plots to determine if there are any visually obvious differences between the two grid densities. There is no real point in doing any quantitative comparisons if pressure and velocity contours are clearly different. Grid refinement should be carried out to the point that visual inspection can no longer tell the difference between contour patterns. Once this has been achieved, actual data must be examined. For the following grid analysis, three levels of refinement are considered. The mesh properties are given in the following table. It should be noted that the portion of the mesh used to resolve the boundary layer does not change with different meshes.

Before looking at derived parameters such as collision potential, it is important to examine variables common among all flows. Pressure coefficient is one such parameter, and it is defined as:

FLUENT allows the pressure coefficient to be plotted, but the user must be careful to set the correct reference values. In this case, p0 = 0 (gauge pressure), ρ = 1000 kg/m3, and V = 0.1 m/s. Figure 2 shows the contours of pressure coefficient for flow around fifteen baffles.

Figure 2: Contours of Pressure Coefficient for 15 Baffles (h/s=10)

Using FLUENT's ability to draw lines/rakes, the pressure coefficient can be extracted and plotted at a constant height. Figure 3 shows the pressure coefficient at a line height of 0.5 m (half the flocculator height for h/s=10) across all the baffle spaces.

Figure 3: Pressure Coefficient at Half the Flocculator Height for all Baffle Spaces (h/s=10)

Finally, the extracted pressure coefficient can be plotted across all baffle spaces for the three different mesh densities (figure 4). A couple of interesting trends can be observed. First, doubling the mesh density (i.e. halving the cell dimension) suggests asymptotic convergence since the difference going from ∆ to ∆/2 is larger than that from ∆/2 to ∆/4. Also, as the flow progresses through the baffles, the difference between the pressure coefficient for each mesh density gets smaller.

Figure 4: Pressure Coefficient Across Each Baffle Space for 3 Different Mesh Densities

To further investigate these phenomena, an average pressure coefficient was calculated for each baffle space for all three meshes, and this is plotted in figure 5.

Figure 5: Average Pressure Coefficient for each Baffle Space for 3 Different Mesh Densities

The average values were used to determine percent differences between Cp for each baffle space with the results appearing in figure 6. This helps confirm the observations mentioned above. As the plot shows, refinement from ∆ to ∆/2 has percent differences ranging from 2% to 2.5% for each baffle space while going from ∆/2 to ∆/4 holds the percent differences below 1%. The percent difference also decreases as the flow goes through each baffle and drops off considerably in the last few spaces. While it is interesting that Cp far from the entrance seems to have decreasing dependence on the cell size, the most important result is probably the fact that the most refined mesh is consistently less than 1% different than the previous refinement. Since the error introduced by the turbulence model usually ranges from 5-10%, it is not necessary to pursue any further increase in mesh density. Thus the ∆/4 grid will be the standard used for all further flocculator investigations.

Figure 6: Percent Difference for Pressure Coefficient in each Baffle Space when Increasing Mesh Density

A quick note on work done in previous years: After reviewing reports going back to the spring of 2008, it appears that the only quantitative look at mesh sensitivity was during the CFD team's first semester of existence. The analysis looked at the pressure coefficient drop around one baffle for three grids that had only a small variance in cell dimension. Hopefully this more systematic grid refinement combined with many baffles and an error analysis will help put to rest any further concerns about grid size.

Iterative Convergence Verification

Iterative convergence is the method used to deal with the nonlinearity of the governing Navier-Stokes equations. Basically the procedure involves beginning with a guess value, then solving the discrete governing equations and comparing that solution to the guess value and finding that difference. The solution then becomes the guess value for the next iteration. When the difference between the guess value and the solution falls below a specified residual, convergence is reached. The error associated with this procedure is based on the residual level. Not only does the iterative process reduce linearization error, it also reduces the error associated with a more efficient (and thus faster) matrix inversion process that is common among CFD codes. Further, residual values can be monitored during a simulation to determine if the solution is smoothly converging or if it is diverging or wildly fluctuating. Figure 7 is a standard residual plot for the five variables of the k-ε realizable model where each variable must converge to 10^-6.

Figure 7: Residuals for the 5 Variables of the k-ε Realizable Turbulence Model

To see the effect of this error on the overall simulation, the pressure coefficient will be examined again. In the following two cases the iterative process was allowed to run until the variables converged to 10^-6 and 10^-9. Figure 8 shows the pressure coefficient plotted across each baffle space. As can be seen, the two levels of convergence are virtually indistinguishable.

Figure 8: Pressure Coefficient Across Each Baffle Space for 2 Different Convergence Residuals

The average values in figure 9 naturally indicate the same result.

Figure 9: Average Pressure Coefficient for each Baffle Space for 2 Different Convergence Residuals

Looking at the percent differences between the two residuals (figure 10) confirms that iterative convergence has a very small contribution to the overall error. While grid error was approximately 1%, iterative error is less than 0.1% even towards the end of the flocculator where results begin to diverge. The bottom line is that 10-6 is more than adequate for iterative residual value.

Figure 10: Percent Difference for Pressure Coefficient in each Baffle Space when Decreasing the Residual

Wall Treatment Verification

Wall treatment is different than other sources of error in that it fits both the verification and validation categories. The verification aspect stems from the grid refinement near a no-slip wall to properly resolve the boundary layer while the validation relates to the specific model used for near wall flows. Here, the discussion will focus on the grid aspect of wall treatment. When looking at a fully developed turbulent flow in a channel, the velocity profile is noticeably different from its laminar counterpart. While laminar flow will have a smooth parabolic velocity profile, the turbulent case will have a mean velocity profile that is much more flattened with a higher velocity gradient near the walls. A finer mesh near the wall will help properly resolve these gradients. Zooming in to focus just on the near wall region, one would find three distinct layers that have different characteristics and mathematical descriptions. A dimensionless measure of distance from the wall (y+) has been developed and the aforementioned layers are bounded by values of this parameter. The three layers of flow near a wall are the viscous sublayer (y+ < 5), the buffer layer (5 < y+ < 30), and the log-law layer (y+ > 30). Beyond these three layers is the core flow and wall treatment no longer applies. The viscous sublayer is dominated by the laminar regime while the log-law layer is almost exclusively turbulent. The buffer layer is where both flow regimes are important, and CFD codes such as FLUENT have difficulty obtaining accurate solutions when the first cell from the wall is in this region. Thus a grid must either be small enough to resolve the viscous sublayer or big enough that it starts in the log-law layer. Since shear gradients from the boundary layer contribute to the turbulent energy dissipation, it is important to have full resolution of near wall flow, so y+ must be kept below five. In previous years, the CFD developed a boundary layer mesh grading that kept y+ small enough, so this year it will suffice to confirm this result. Figure 11 is wall y+ for a boundary layer grid that is nine cells with the first cell 0.0003 m thick and growing at a rate 1.25 per cell. As can be seen, the value of y+ is completely below 5 for all locations along the wall indicating the viscous sublayer is fully resolved, so there is no extra error associated with having the first cell in the buffer layer.

Figure 11: Wall y+ throughout the Flocculator

Validation

Even after minimizing the error from the mesh, the wall mesh, and iterations, the turbulence model may still not be accurately modeling flow in the flocculator and, more generally, around a 180° bend. When the project first began in 2008, the k-ε standard model was originally decided on because it was the best match to the observed flow in the demonstration plant. Following the realization that the demo plant was in the laminar regime, a more scientific (though probably still incorrect) analysis was undertaken in the fall of that year to select a turbulence model. The well-studied flow over a back step was used as a benchmark to determine that the k-ε realizable model was actually better for modeling the recirculation and reattachment around each baffle. Presently this is the model used for flocculator simulations. Unfortunately, the back step does not account for the turning flow behavior prior to the separation point, so it may not be fully adequate for determining the best model. The primary reason for skepticism of the turbulence model comes from experimental data obtained from a test flocculator. A team member from a separate subteam built a model flocculator with a height of 0.3 m, width of 0.6 m, and spacing and clearance height of 0.07 m. The inlet velocity was 0.103 m/s, and pressure drop data was recorded. The experimental results yielded a value for the pressure coefficient drop (K) of 2.29. This value is contrasted with a K value of ~3.9 for the same geometry and inlet conditions from a CFD simulation. Since grid, iterative, and wall treatment error have all previously been addressed, it appears only the turbulence model can account for this variation. Based on documentation of experimental data from FLUENT showing the Spalart-Allmaras and Reynolds Stress models could best handle a 180° turn (figure 12), these two models were tested with the same settings as the flocculator experimental setup with the hope of finding a replacement turbulence model.

Figure 12: Normalized Velocity for different Turbulence Models around a 180° Turn (FLUENT, 2001)

The Spalart-Allmaras was tested first since it is a one equation turbulence model and thus less computationally expensive. While this model did converge quickly, clearly the k-ε realizable is still the better model since it is much closer to the experimental value for all baffles. This can be seen in figure 13 where K is plotted in each baffle space for both models. This simple check helped to eliminate the one equation model from consideration.

Figure 13: Kbaffle Comparison of k-ε realizable and Spalart-Allmaras

The other possibility was the Reynolds Stress model which appeared to have the greatest potential for successfully modeling the 180° turn based on the data from FLUENT, but it is a multi-equation model and takes longer to run. It was also discovered that this model is very sensitive to the initial guess since it had trouble converging during the initial attempts to use it for modeling the flocculator. In order to achieve convergence, the simulation had to initially be run with the k-ε realizable model until it converged to 10-6, and then the model was switched over to the first order accurate Reynolds Stress model which was also run to 10-6. Attempts to then run the second order accurate solver were unsuccessful and will need to be looked at next semester, but from the first order solution, it appears that the Reynolds Stress model does indeed do a better job of predicting pressure coefficient drop (figure 14). Clearly validation work will need to continue into next semester as the CFD and experimental values are still way too far off to be acceptable.

Figure 14: Kbaffle Comparison of k-ε realizable and Reynolds Stress

Flow Field Description

Following the lengthy verification and validation procedure, simulations can be run to obtain information about the flow properties and performance parameters. Since AguaClara is looking for a complete understanding of flow through the flocculator, it is the CFD team's task to come up with the model for this purpose. The entire flocculator clearly cannot be modeled, meshed, and solved as this would take a prohibitively large amount of time, so some portion of the flocculator will need to be modeled that is representative of most of the flow through the baffles. It was originally thought that there were two ways to do this, and one method involved using a periodic solver that would only model one baffle turn and impose the requirement that flow variables be the same at the inlet and outlet. While this method was ultimately abandoned, the following is a quick overview of this model and an explanation for why it did not work.

Periodic Case

The periodic simulation was a proposed method of determining flow characteristics far downstream of the inlet. Instead of creating geometry for a large number of baffles and looking at the solution near the end of that case, only one baffle would need to be created with the boundary conditions providing the necessary information to the solver. The time it takes to create and mesh a multi-baffle flocculator can be extensive, so clearly being able to examine this case using only one turn would save a significant amount of work. Time would also be saved during the solution step of the procedure. For each iteration, FLUENT solves the discrete governing equations in each cell of the mesh. The more baffles a particular case has, the more cells are required to obtain an accurate solution. A large number of cells can significantly slow the time each iteration takes, even with the powerful processors available in the Swanson lab. While a one baffle turn might take about a second per iteration, ten baffles might take ten seconds per iteration, and since it takes thousands of iterations to reach the desired residual level, these simulations could end up taking more than a day to compute. For the above mentioned reasons, a periodic baffle is desirable, and much effort was placed on implementing this case.

According to the FLUENT user guide, the following are conditions and limitations on its periodic solver:

• Flow must be incompressible
• The geometry must be translationally periodic
• No extra mass sources/sinks
• No reacting flow
• No multiphase flow

The periodic baffle case meets all these criteria, so a solution should be obtainable.

Boundary conditions for the periodic case can only come in two forms. Either a prescribed mass flow rate or pressure gradient. The mass flow rate was attempted first, and a value of 6 kg/s was used corresponding to some of the lower flow rate plants (it should have been 10 kg/s for this particular case but it didn't matter since no solution could be found). After several attempts with different under-relaxation factors and thousands of iterations, the general pattern was a slow and steady convergence of the two components of momentum and the two turbulent parameters contrasted by a sharp divergence in continuity. In other words, the solver was not finding a solution. A shift to specifying a pressure gradient of 15 Pa/m did little to improve this situation. The five variable residuals would all approach residual levels of 10^-4 then begin to fluctuate as the solver dealt with instabilities in values between iterations. Ultimately it was determined that the flow had to be periodic along the axis of translation which is not explicitly stated in the guide. This was confirmed after communicating with FLUENT technical support. As such, the periodic case was abandoned around the middle of the semester in favor of modeling for baffles.

More Baffles

The failure of the periodic case created the motivation for modeling more baffles. While it takes more time to draw the geometry, create the mesh, and run the simulation, it is necessary for fully describing flow through the flocculator. By modeling more baffles, it is possible to determine if and when the flow around each turn becomes "baffle-wise" periodic. This simply means that after the flow has gone around enough baffles, it becomes independent of entrance effects and starts repeating in each baffle space. Figure 15 provides a qualitative look at what is meant by baffle-wise periodic. As can be seen, the first baffle space has very little turbulent energy dissipation, the second baffle space has a larger area of turbulent energy dissipation, and the third baffle space and beyond have similar looking energy dissipation zones. This very crude visual inspection would indicate that the flow is baffle-wise periodic after only two turns. Obviously, visual inspection is not sufficient to determine when flow in each baffle space starts repeating, and a more rigorous quantitative analysis would need to be used. Still, it provides a useful definition of the concept.

Figure 15: Contours of Turbulent Energy Dissipation, ε, (m2/s3) through 15 Baffles (h/s=10)

Once this repeating flow has been achieved, it can easily be extrapolated through the entire flocculator providing the full flow description that AguaClara requires. Prior to this year, more than five baffles had never been modeled. While this may have satisfied the qualitative procedure mentioned above, for quantitative purposes, the number of baffles needed to increase. The current mesh density combined with the limit on the number of cells allowed with the FLUENT educational license capped CFD simulations to approximately fifteen baffles, so this became the new standard for investigations. Figures 16 and 17 are contours of turbulent energy dissipation for smaller h/s ratios. It is important to look at varying h/s values because the number of baffles needed to reach baffle-wise periodic flow is different.

Figure 16: Contours of Turbulent Energy Dissipation, ε, (m2/s3) through 15 Baffles (h/s=2)

Figure 17: Contours of Turbulent Energy Dissipation, ε, (m2/s3) through 15 Baffles (h/s=4)

User defined functions can be used to extract performance parameters of interest from the simulation data, and the primary performance parameter is collision potential. The θ variable in the equation is cell area (assuming unit depth for a 2D simulation) divided by the volumetric flow rate. Collision potential was extracted and plotted for all baffles spaces for h/s values of 2, 4, and 10. The results are presented in figures 18, 19, and 20 respectively.

Figure 18: Collision Potential v. Baffle Space # (h/s=2)

Figure 19: Collision Potential v. Baffle Space # (h/s=4)

Figure 20: Collision Potential v. Baffle Space # (h/s=10)

In order to get a more refined look at Ψ, the first few baffles for each case will be ignored since they are clearly dominated by entrance effects. Eliminating the first four baffle spaces for each case yields the plots shown in figures 21, 22, and 23.

Figure 21: Collision Potential v. Baffle Space #'s 5-16 (h/s=2)

Figure 22: Collision Potential v. Baffle Space #'s 5-16 (h/s=4)

Figure 23: Collision Potential v. Baffle Space #'s 5-16 (h/s=10)

The three plots above have the same y-axis range, so it appears the h/s of 10 reaches baffle-wise periodic flow after fewer baffles than 4 and 2. Removing four more baffles will lead to even further refinement, and these plots are given in figures 24, 25, and 26.

Figure 24: Collision Potential v. Baffle Space #'s 9-16 (h/s=2)

Figure 25: Collision Potential v. Baffle Space #'s 9-16 (h/s=4)

Figure 26: Collision Potential v. Baffle Space #'s 9-16 (h/s=10)

It appears that after about nine baffle spaces (eight baffles) the flows for the h/s values of 4 and 10 have reached the desired periodic behavior. The smaller h/s of 2 still appears to need more baffles to determine when this type of flow occurs, but for the majority of flocculator investigations where h/s is greater than 2, fifteen baffles is more than adequate to for determining values of collision potential far from the entrance.
One quick note about the plots above: In the last baffle space for all three h/s ratios, the collision potential differs significantly from the periodic value in the preceding baffle spaces. For h/s values of 4 and 10, Ψ drops in the last baffle space because there is no turn at the flocculator exit, so the only main contribution to energy dissipation comes from the flow around the final baffle. The reason it increases for h/s=2 has to do with flow reattachment not occurring before the flow reaches the exit. FLUENT has difficulty modeling flow in this situation, and this difficulty likely causes the much higher ε value. In all cases, the final baffle space data will be ignored since it is not representative of flow through the majority of the flocculator.

Conclusion

Overall, after a month or so spent learning the FLUENT-Workbench interface and working out issues related to geometry sketching and meshing, much of the semester consisted of the tedious but necessary verification and validation process to help ensure the results obtained from flocculator simulations are correct. Grid error has been effectively reduced to less than 1%, and iterative error was found to be insignificant at less than 0.1%. Further, the near-wall mesh is fully adequate for resolving the viscous sublayer since it keeps y+ values below five. The k-ε realizable turbulence model has been shown to not accurately predict the drop in pressure coefficient around each baffle, and the Reynolds Stress model seems like the best candidate to serve as a replacement. To obtain a flow field solution through the entire flocculator, it was necessary to increase the number of baffles to fifteen since the periodic case did not work and five baffles was not adequate for ensuring repeating flow through each baffle had been achieved. For the end of this semester, the goal will be to complete the turbulence model validation process. Once this is done, accurate performance parameters can be extracted and used as plant design considerations.