Computational Fluid Dynamics

Flocculation Tank Simulation

Overview

The FLUENT analysis software used in the Fall 2008 and Spring 2009 semesters was changed over the summer of 2009. FLUENT was integrated into the ANSYS software package. This changed the procedure for constructing CFD simulations. The Gambit interface for building meshes was removed and mesh construction is now handled in Workbench. The sketching of geometries is handled in a SolidWorks style environment and the sketches are transferred into a meshing tool. Analysis of fluid flow is a still handled by the classic FLUENT solver.

The shift from Gambit to ANSYS for geometry sketching rendered the geometries and automated mesh generating journal files in Gambit, from Fall 2008 and Spring 2009, incompatible with ANSYS, unless the files were converted with a MATLAB code. New 2D sketches of AguaClara baffle systems were created in ANSYS. Creation of lines from sketches and the line projection command are essential for creating systems of baffles in the new version of Fluent. For turbulent flow in FLUENT, the mesh must be refined enough for the solution of the simulation to stay out of the mixed boundary layer regime of turbulent y+ = 5 to 30. The turbulent y+ in the solution can be below 5 or above 30 if the cells at the walls of the mesh are made small enough to solve the turbulent boundary layer gradients near the wall. This is achieved with the biasing tool in the mesh creation interface. The mesh can be biased to have smaller cells near the wall than at the center of the pipe for example.

Unable to render embedded object: File (mesh biasing.png) not found.

Biasing of mesh near the wall

The investigation began with setting up simulations from Fall 2008 in the new FLUENT interface. We ran these cases to verify the new version of FLUENT. The most extensive investigation was in periodic boundary conditions. This tool in the FLUENT solver iterates the simulation until convergence with the requirement that inlet and outlet values are the same. This tool is capable of performing this analysis strictly for symmetric geometries. A two turn baffle model was used for the periodic analysis. Periodic solution allows the user to initialize the flow at the inlet with the mass flow rate (kg/s) or pressure drop/length traveled (Pa/m). Cases of 5, 10, and 15 baffles were modeled and simulated to compare with the converged periodic solution. The user defined functions (UDFs) were used to explore how to analysis the performance parameters in the said cases.

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.

FLUENT Updates

Recently ANSYS acquired FLUENT and integrated it with the more user-friendly Workbench software. While this makes new projects much easier to set up when compared to the previous Gambit interface, it required a good deal of work to adapt previous work to the new system. Workbench is not a new tool itself, and ANSYS has used it for a number of years for its structural finite element solver. It allows the user to step through the process of geometry creation and meshing, setting up the physics, solving the problem, and postprocessing the results. Work began immediately on reproducing results from past years with the goal of ensuring the new system did not fundamentally change any of the results that had previously been obtained. It was clear early on that, while better than building geometry and meshes in GAMBIT, the new system left something to be desired. Certain common sense tasks that someone would likely want to perform in any CAD style program required nonintuitive methods. For example, when creating a baffle turn in the flocculator, our 2D model was best served by having an infinitely thin wall where the fluid would flow in opposite directions on different sides. The software help documentation was useless for providing an answer, and only after consulting technical support did we find that lines had to be drawn, created as separate bodies, and then projected onto the surface through which fluid would flow. Another such issue occurred while creating the near-wall mesh. In order to properly resolve the viscous sublayer (i.e. y+ values less than 5), cell size growth had to occur from an initial cell size (0.0003m) to the size used for meshing the flow's core that is "far" from the wall. This is implemented using a bias factor that governs the change in size of each cell when moving away from the wall. Unfortunately, there is no simple user input where first cell size, number of divisions, and total length are specified and the bias factor is simply computed from this information (Gambit actually had this ability which indicates ANSYS took this a step backwards in my opinion). Also the documentation doesn't provide this algorithm leaving it to the user to either figure it out independently or simply guess until the correct size is obtained.

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 pressure coefficient at a height of 0.5 m (half the flocculator height for h/s=10) across all the baffle spaces for the three mesh densities. 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 2: 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 3.

Figure 3: 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 4. 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 4: 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 5 is a standard residual plot for the five variables of the k-ε realizable model where each variable must converge to 10^-6.

Figure 5: 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 6 shows the pressure coefficient plotted across each baffle space. As can be seen, the two levels of convergence are virtually indistinguishable.

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

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

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

Looking at the percent differences between the two residuals (figure 8) 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 8: 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. Wall y+ is given as:

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

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.

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, but it now appears this method will not work.

More Baffles

The motivation for more baffles came from the inability of the periodic solution to converge. To advance the project, it was essential to find the flow solution after enough baffles to achieve time independent conditions. While the simulations would take much longer, the decision was made to proceed with more baffles and look at the relative change in flow variables before and after successive turns. The solution took ~3000 iterations (6+ hours) to converge, but the process was worth the time since the results seem to indicate that after only about 5 to 10 baffles, the solution becomes steady.

Contours of Velocity (m/s) over 15 Baffles (H/S = 10)

Looking at turblent energy dissipation also revealed the same sort of constant values after 5 or so baffles.

Contours of Turbulent Energy Dissipation (m^2/s^3) over 15 Baffles

Conclusion

Multiple trials were tested to attain a periodic solution. The mass flow rate input was calculated to be 10 kg/s, based on a density of 1000 kg/m^3, an area of 0.1 m X 1 m, and a velocity of 0.1 m/s. Professor Weber-Shirk estimated a pressure drop of 30 Pa over two baffle turns, which was used in the pressure drop initial condition for periodic solution. The Fluent solver solution controls were adjusted to allow the periodic simulation to converge. The behavior of the converged periodic solutions could not be validated. The mass flow rate initialized periodic solution had a very large pressure drop over the baffle turns. The pressure drop periodic boundary conditions yielded little if no pressure difference between the outlet and the inlet. To this point the periodic solution for the AguaClara baffles does not work. The 10 and 15 turn baffle solutions show stability of energy dissipation rate after multiple turns, perhaps showing promise for analyzing the performance parameters for optimal baffle geometries far down stream. The UDFs are being adapted to analyze the current cases.

Post-Midterm Analysis

Flocculator Design Case