Towards a Dynamic Compartmental Model of a Lamellar Settler

Abstract

Lamellar settlers are widely used in solid separation due to their relative high surface area paired with a small tank volume. Settlers exposed to high load variations may change their settling properties based on the variation of the flow pattern. Applying the critical point theory and coherent structure concept, a flow topology analysis was performed with the help of computational fluid dynamic simulations. The compartments were determined by the fluid flow, and the dynamic behavior of the compartments was taken into account. Under normal diurnal load variation, the architecture of the compartments did not change, in contrast to the mass transport between the zones, whereas the sludge removal process made significant changes in the architecture. The results showed two main flow zones with significant internal recirculation in the first zone. The hydraulic surface loading was examined in each zone, and the study revealed that the actual hydraulic load could be from 1.5 to 4.5 times much higher than that in individual zones due to recirculation eddies in the reactor tanks. The design hydraulic loading rate did not consider the local flow pattern. The discrete phase model approximation provided acceptable results, and the extent of the recirculation zone changed stepwise with different loads.

1. Introduction

The removal of solid matter, including both organic and inorganic components, is a crucial aspect of water and wastewater treatment processes. This solid separation step is responsible for reducing the pollutant load of downstream equipment by removing a significant portion of the contaminants present in raw wastewater [1]. The efficiency of this separation process is dependent on various factors, including the characteristics of the raw water, such as the solid content and particle size distribution, as well as the hydrodynamic conditions and environmental factors present [2].

From a hydrodynamic perspective, the overall efficiency of the separation process is influenced by several key factors, including the hydraulic load, the density difference between the solid and the liquid phases, and the viscosity of the liquid. Additionally, the gravitational acceleration is also an important factor to consider. At the local level, the presence of turbulence and shear forces could also play a role in the efficiency of the solid–liquid separation process [3]. The shear rate required for optimal separation can vary depending on the characteristics of the particles being separated. Flocculating particles, which tend to aggregate, require a lower shear rate to allow for particle aggregation and settling than non-flocculating particles, which do not aggregate as easily. Additionally, the desired shear rate can vary depending on the specific stage of the separation process: a low shear rate is desirable during the sedimentation stage to allow for the particles to settle out of the liquid. Undesirable shear can negatively impact the efficiency of the separation process by causing the formation of small flocs which are more difficult to settle, destabilizing emulsions and decreasing the efficiency of flocculating agents. Additionally, high shear rates can cause the particles to break apart, leading to a decrease in separation efficiency [4].

Local turbulence parameters such as turbulent viscosity and turbulent dissipation have also an effect on the phase separation efficiency. Turbulent viscosity originates from the rate of strain in the flow; with a high virtual viscosity, the flow resistance to the settling motion increases, and thus the separation efficiency decreases. In contrast, high turbulent dissipation refers to a high local mixing intensity and promotes flocculation and settling [5].

Furthermore, the flow field topology within a process unit could provide information on the hydrodynamic behavior. A topology that promotes the formation of eddies and vortices could increase the particle collision and, thus, the flocculation, which may lead to higher settling rates [6]. An opposite effect could be observed if the flow is highly turbulent, or shear layers are developed [7]. The flow topology can change over time as the flow evolves and the conditions in the liquid change [8]. The fluid structure within a sedimentation process can be affected by variations in the load conditions, such as changes in the flow rate, particle size, and concentration, and properties of the liquid. The optimal design of the sedimentation process should take into account these potential variations and adjust accordingly. Such adjustments may include the geometry of the tank, with a focus on the radial flow, the positioning of inlets and outlets, the implementation of baffles, and the incorporation of inclined plates.

The use of inclined plates increases the total surface area available for solid–liquid separation, thus reducing the required size of the process unit while maintaining efficiency [9]. The removal of settled sludge can be performed in a continuous or in a batch manner. The discharge flow from the sludge removal process may also have an impact on the flow pattern through backpropagation [10].

Computational fluid dynamics (CFD) could be a helpful tool to reveal the possible small-scale flow phenomena within the process unit and optimize the hydrodynamic conditions with variable loads. Many efforts have been made in the numerical modelling of the fluid flow in lamella settlers [11]. Earlier studies widely used single-phase flow models because they were computationally efficient and produced accurate results [12]. With recent advancements in computing power, researchers have shifted to using multiphase models to more accurately model sedimentation [13]. The Euler–Euler or Euler–Lagrange method can be used to simulate the multiphase flow observed in settling tanks. In the near future, a full-scale, realistic transient simulation of a settling tank that includes sub-models for factors such as flocculation, sludge rheology, and thermal conditions will be possible [14].

As for the transport processes, the approaches to describing reactors range from the simplest, treating the entire reactor as an ideal one, to the most detailed, using CFD with fine discretization. The intermediate approach considers the reactor as a combination of some ideal reactors, the so-called systemic approach or the compartmental method, where all compartments are assumed to be perfectly mixed. The compartments can be defined purely by geometry or by a certain flow characteristic [15].

Coherent structures in a flow are patterns that are distinguished from the surrounding turbulent or chaotic flow by their relative stability and regularity. These can be characterized by strong gradients in velocity and/or vorticity [16]. Coherent structures in a flow can exhibit both symmetry and asymmetry. It depends on the specific flow and the structures that are observed within it. Some coherent structures, such as vortices or jets, can be symmetric in shape or behavior, while others can be asymmetric. For example, a symmetric vortex will have a circular shape, and the fluid motion within the vortex will be rotationally symmetric about its central axis, while an asymmetric vortex will have an irregular shape, and the fluid motion within it will not be rotationally symmetric. It is also possible for a flow to have both symmetric and asymmetric coherent structures present at the same time [17]. The identification could apply the linear combination of the proper orthogonal decomposition (POD) modes of vorticity [18] or criterion-based assessments [19].

Based on the coherent structures detected in a flow field, zones with similarity could be delineated. There are several techniques to separate compartments. The most common is the visualization approach, where the streamlines delimit the regions, whereas a spectral analysis decomposes the flow into various frequency components and identifies the coherent structures that are associated with a certain frequency range [20]. While Lagrangian methods take into account particle tracking, the Eulerian approach applying a topological analysis originates from fixed-point theories [21]. Both are widely used in steady-state conditions, but the dynamic behavior of the delineated zones is rarely investigated. There are examples of Artificial Intelligence (AI)-based methods in vortex detection [22], especially originating from their image recognition capability, but the mechanistic models have an advantage in terms of robustness and transparency.

The objective of this study was to predict the hydrodynamic conditions of an inclined/lamellar settler utilizing compartmental modeling and to study the compartments’ configuration and inter-compartmental mass transfer at various hydraulic loads. Although CFD-based research is not novel in this field, the innovation of this study lies in the compartmentalization approach and the consideration of the compartments’ dynamics. The compartments were defined using critical point theory and dynamic behavior, and the evolution of the volume of the compartments based on the hydraulic load incoming to the settler was investigated, which to the authors’ knowledge, has not been investigated to date.

2. Materials and Methods

2.1. Flow Field Calculation

The governing equations of CFD simulations are usually based on the Navier–Stokes equations, which describe the motion and behavior of fluids. These equations can be expressed as a set of partial differential equations (PDEs) and can be solved using numerical methods. For 3D steady-state incompressible flows, applying Reynolds averaging and Einstein notation in Cartesian coordinates, the following form was used for continuity Equation (1) and the three components of the momentum equations (Equation (2)).

$$\frac{\partial {\overline{u}}_i}{\partial x_i}=0$$

(1)

$$\rho \overline{u_j}\frac{\partial {\overline{u}}_i}{\partial x_j}=\rho {\overline{f}}_i+\frac{\partial }{\partial x_j}\left. [-\overline{p}{\delta }_{ij}+\mu \left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}\right)-\rho \overline{u_i^{\prime }{u}_j^{\prime }}\right]$$

(2)

where i, j: x, y, z are the Cartesian coordinates, u is the velocity, ρ is the fluid density, f is the external force, g is the gravitational acceleration, p is the pressure, δ is the Kronecker delta, µ is the dynamic viscosity.

The left-hand side of Equation (2) contains the convective term, while the right-hand side contains the external force, pressure gradient, diffusion term, and virtual stress from the Reynolds averaging. For turbulence closure, the k-ω model was applied, which assumes isotropy and near-wall treatment and is well-suited for flows with strong pressure gradients, such as those found in boundary layers and separated flows [23]. As for the numerical approach, the discretization scheme followed the upwind scheme. Second-order Taylor-series expansion was applied in order to reduce the numerical diffusion [24]. The general SIMPLE scheme (Semi-Implicit Method for Pressure-Linked Equations) was solved for pressure–velocity coupling. Steady-state flow field separate simulations for various incoming flows were performed. The convergence criteria were to achieve that the residuals of the governing equations between two iterative steps were minimized.

Lagrangian particle tracking in the Eulerian flow field involved following the movement of particles in a fluid as they were carried along by the flow. This was carried out by tracking the position and velocity of each individual particle over time. There was only a one-way coupling between the solid particles and the fluid flow; the particles were assumed to be small enough that their motion could be approximated as being affected by the surrounding fluid, but not significantly altering the fluid’s flow. The random walk approach was used to capture the turbulent properties [25]. The mathematical formulation of Lagrangian particle tracking follows the force balance on the particle:

$$m_p\frac{d{u}_p}{dt}=\frac{1}{2}{\rho }_f{C}_D\left|\overline{u}-u_p\right|\left(\overline{u}-u_p\right)A_p+\left({\rho }_p-{\rho }_f\right)g{V}_p$$

(3)

where m_p is the particle mass, V_p is the particle volume, A_p is the particle surface perpendicular to the movement, u_p is the particle velocity, u is the flow velocity, ρ_p is the particle density, ρ_f is the fluid density, C_D is the drag coefficient calculated as follows:

$$C_D=\left\{\begin{array}{c}\frac{24}{Re}\left(1+0.15R{e}^{0.687}\right), Re<1000\\ 0.44, Re>1000\end{array}\right.$$

(4)

where Reynolds number (Re) can be determined as follows:

$$Re=\frac{{\rho }_f{d}_p\left|\overline{u}-u_p\right|}{\mu }$$

(5)

where dp is the particle diameter.

To predict the dispersion due to turbulence, the Gaussian-distributed random velocity fluctuations were calculated as follows:

$$u_i^{\prime }=\vartheta \sqrt{\overline{u_i^{\prime }{}^2}}$$

(6)

where $\vartheta$ is a normally distributed random number. The characteristic eddy lifetime (τ_e) was then calculated applying ω from the turbulence model, as follows [25]:

$${\tau }_e=\frac{0.3}{\omega }$$

(7)

The methodology of the residence time distribution (RTD) analysis was based on the author’s previous work [26], in which a two-phase flow was assumed, where both phases of the fluid had the same material properties, and the only difference between the phases was their spatial and temporal distribution. As an initial condition, the entire domain was full with phase 1 using an a priori steady-state flow simulation. At a given time (t = 0 s), the simulation was set to transient, and the input boundary condition was changed so that only phase 2 was discharged at the same rate as phase 1. As time progressed, phase 1 was displaced because there was no refill through the input boundary condition, and the rate of this process could be tracked by the volume fraction (α) scalar variable:

$$V_q=\iiint {\alpha }_{q}dV$$

(8)

$${\sum }_{q=1}^n{\alpha }_q=1$$

(9)

where α is the volume fraction of the q-th phase.

The volume fraction as a function of time could be depicted at the outlet boundary, which represents the so-called F-curve in RTD analysis [27]. The cumulative function of the volume fraction applied the normalized time as follows [28]:

$$\theta =\frac{t}{\overline{t}}$$

(10)

where $\theta$ is the normalized time, and $\overline{t}$ is the average time.

2.2. Coherent Structure Detection and Compartmentalization

The critical point theory provides a useful framework for understanding the complex behavior of fluid flows with coherent structures. The critical point is a point where the gradient of the velocity field is zero. The velocity field changes from a region of acceleration to a region of deceleration, leading to the formation of vortices and other coherent structures in the flow. By the help of the local velocity gradient tensor, the behavior of the fluid in the near zone of critical points could be outlined by applying the eigenvalues and eigenvectors of the matrix. The Jacobian matrix for a two-dimensional vector (u, v) is

$$\mathit{A}={⌈\frac{\partial u_i}{\partial x_i}⌉}_{x_{i,0}}={\left. [\begin{array}{cc}\frac{\partial u}{\partial x}& \frac{\partial u}{\partial y}\\ \frac{\partial v}{\partial x}& \frac{\partial v}{\partial y}\end{array}\right]}_{x_0,y_0}$$

(11)

The eigenvalues of the A matrix can be calculated using the following equation:

$$det \left(\mathit{A}-\lambda \mathit{I}\right)=0$$

(12)

where λ is the eigenvalue, and I is the unity matrix.

If the real part of the eigenvalue is negative (Re < 0), it creates an attractive force, whereas if it is positive, it causes a repulsive force (Re > 0). On the other hand, the imaginary part (Im) of the eigenvalue results in a rotation of the vector field around the critical point. If there is an imaginary eigenvalue, the vector field around the critical point will be entirely rotational. This categorization is based on the existence of the imaginary part and the sign of the real part of the complex eigenvalue based on

Table 1. Classification of critical points based on the eigenvalue of the velocity gradient tensor.

Name of the Critical Point	Real and Imaginary Part of the Eigenvalue
Saddle point (S)	Re1 < 0, Re2 > 0, Im = 0
Attracting node (AN)	Re1 < 0, Re2 < 0, Im = 0
Repelling node (RN)	Re1 > 0, Re2 > 0, Im = 0
Attracting focus (AF)	Re1 < 0, Re2 < 0, Im < >0
Repelling focus (RF)	Re1 > 0, Re2 > 0, Im < >0
Center (C)	Re1 = 0, Re2 = 0, Im < >0

[29].

By examining a steady-state flow field, streamlines were used to identify critical points on visual and intuitive construct compartments. The separatrices served as boundaries for each compartment, and in turbulent flows, mass exchange occurred between zones. The fluxes through the compartments were calculated to establish the architecture of the compartment model, based on a one-phase steady-state flow.

2.3. Mesh Independence Test and Validity of the Model

The goal of the mesh independence analysis was to investigate whether the further refinement of the grid size could significantly improve the accuracy of the results. The test started with a coarse grid and proceeded with a gradual increase of the number of cells or elements in the mesh. As a result, coarse, practical, and fine meshes were created with cell numbers of 87,126, 348,523, 1,349,220, respectively (see details in Figure 1).

The error (ε) could be estimated as follows:

$$\epsilon =\frac{{\phi }_{fine}-{\phi }_{practical}}{r^p-1}$$

(13)

where φ is the characteristic variable for the flow, r refers to the aspect ratio of the mesh element size, and p is the order of the convergence calculated as follows:

$$p={\log }_r\frac{{\phi }_{practical}-{\phi }_{coarse}}{{\phi }_{fine}-{\phi }_{practical}}$$

(14)

2.4. Model Setup

Inclined plates were placed in a rectangular channel in order to increase the settling area. The basis of the geometry and the arrangement was published in Dao et al. (2019) [30], taking into consideration the erratum by the authors published in 2020 [30]. The simulation domain had a length of 4.0 m, a height of 2.0 m, and a width of 1.0 m. Eight lamellae were placed in the simulation domain as Figure 2 shows. As for the boundary conditions, mass flow at the inlet and atmospheric pressure at the outlet were set. The inclined plates were walls, the sides of the channel had symmetry boundaries, and besides the free surface, the rest bordering surfaces were no-slip walls.

The hydraulic surface loading rate (HLR) can be calculated by dividing the flow discharge by the total surface area of the settling tank surface and the horizontal projection area of the lamellae as follows:

$$HLR=\frac{Q}{A_0+n{A}_L}$$

(15)

where Q is the flow discharge, A₀ is the settling tank surface area, n is the number of the lamellae, A_L is the surface area of one lamella. Steady-state runs were calculated for various scenarios based on the hydraulic surface loading rate (see

Table 2. Steady-state scenarios.

Scenario Name	Hydraulic Loading Rate [m/h]
Dry weather diurnal average (DWDA)	1
Dry weather diurnal hourly peak (DWDP)	2
Dry weather diurnal hourly low (DWDL)	0.5
Wet weather instantaneous (WW	4

), taking into account the actual operational conditions of a lamellar settler.

Once the Eulerian flow fields had been calculated for each scenario, particles were introduced into the flow, and their movement was monitored. The settling velocity of the particles ranged from 0.24 to 0.80 mm/s, with an even distribution as described in reference [30]. At the boundaries, the particles either escaped through the outlet, which was reflected in the wall shadow zone, or became trapped in the wall zone of the lamellae or at the bottom [31]. The efficiency of particle separation was determined using the following formula:

$$\eta =\frac{N_{inject}-N_{unfinished}-N_{escape}}{N_{inject}-N_{unfinished}}$$

(16)

where η is the particle separation efficiency, and N is the number of particles. The number of unfinished particles refers to those particles whose fate could not be determined within the designated timeframe. Literature data were used to verify the flow model and compare the efficiency of particle transport and sedimentation.

3. Results and Discussion

3.1. Flow Field and Mass Transport Calculation and Validation

A grid independence analysis was conducted using the average flow during dry weather condition. The integral-averaged velocity magnitudes were determined to be 2.92 × 10⁻⁴, 2.85 × 10⁻⁴, and 2.83 × 10⁻⁴ for the coarse, practical, and fine meshes, respectively. As a result, the error was around 10⁻⁸, which could be considered negligible and is likely attributable to rounding.

For each of the four loading scenarios, an analysis of the residence time was conducted, and the surface average volume fraction at the outlet was measured. As time progressed, an increasing amount of fluid from phase 2 reached the outlet. Figure 3 shows the cumulative function (F-curve) for each scenario. The x-axis displays the normalized time.

Starting from time zero, a tracer introduced into the system appeared at different times at the outlet section depending on the hydraulic surface load. For a hydraulic loading rate (HLR) of 1 m/h or greater, it represented 0.36–0.4 of the average residence time, while at low flow rates, it reached a value of approximately 0.54. This was due to slower water movement, lower turbulence, and more dominant plug flow. At the outlet section, the original fluid and tracer could be measured in roughly equal proportions after about 50% of the average residence time had elapsed (HLR = 2 m/h) and in 65% of the other loading scenarios. After the average residence time had passed, the detectable tracer phase ratio varied widely, ranging from 0.68 to 0.83, indicating that the flow pattern varied depending on the HLR. At twice the average residence time, the tracer ratio was between 0.88 and 0.92, while at relative residence times greater than 2.5, the detectability of phase 2 was similar regardless of the HLR. The formation of stable dead zones might account for this, which could not be eliminated even with higher water flow rates. Mass and energy transport between the dead zones and the main flow were solely due to turbulent diffusivity. The study conducted by Dao et al. [30] found that the tracer reached the outlet at θ = 0.35, while the value at θ = 0.53 was associated with alpha = 0.1, indicating good consistency between the present study and the reference cited.

Lagrangian particle tracking was applied to the a priori calculated Eulerian flow field to estimate the separation efficiency. In other words, discrete phase modelling was applied as detailed in the Materials and Methods section. The separation efficiency was compared with the work of Dao et al. [30], who had used the algebraic slip multiphase model and treated the solid phase as a continuum. As shown in

Table 3. Particle separation efficiency.

	Actual Simulation	Literature Data [30]	Deviation
DWDL	0.94	no data
DWDA	0.84	0.88	5%
DWDP	0.6	0.65	8%
WW	0.56	no data

, where the comparison was available, the two different approaches were in reasonable agreement. The deviation was calculated from the difference between the literature data and the actual simulation results, normalized by the actual simulation.

Lagrangian particle tracking underestimated the separation efficiency, which might be due to the boundary condition set at the lamellar zone, since the shadow zones could reflect the total particle mass back to the fluid flow.

3.2. Coherent Structure Detection

To detect the coherent structure, the first step was to reduce the dimensionality of the flow. This involved comparing the volume-averaged velocity magnitudes and the velocity component perpendicular to the main stream in each scenario, as shown in

Table 4. Average velocities and turbulence-related parameters at the simulation domain at various hydraulic loads.

	DWDL	DWDA	DWDP	WW
velocity magnitude (m/s)	1.16 × 10−4	2.85 × 10−4	6.9 × 10−4	1.21 × 10−3
perpendicular velocity (m/s)	6.52 × 10−8	9.47 × 10−8	1.38 × 10−7	1.88 × 10−7
turbulent kinetic energy (m2/s2)	3.33 × 10−10	8.7 × 10−10	5.58 × 10−9	7.46 × 10−8
turbulent dissipation rate (m2/s3)	9.43 × 10−13	1.52 × 10−12	7.5 × 10−12	1.53 × 10−11
turbulent viscosity (kg/ms)	3.86 × 10−6	1.18 × 10−5	7.03 × 10−5	6.4 × 10−4

.

As shown in Table 4, the velocity magnitude was nearly four orders of magnitude higher than the velocity component perpendicular to the flow direction. Cross-stream dispersion was negligible and did not increase linearly with the increasing load. This implies that, for creating compartment models, it was sufficient to project the streamlines onto a two-dimensional plane and determine the critical points based on that. The turbulent kinetic energy and the dissipation rate increased with the hydraulic loading rate, together with the higher turbulent viscosity, leading to a higher dispersion of the particles.

The goal was not to map out all critical points available, but to locate the locations that were significant for settling and material flows. Therefore, the steady-state flow pattern was examined with the fluid entering the domain. Critical points were identified based on the visualization of the streamlines. A critical point occurred directly at the inflow, indicating a small vortex zone, and there were stagnation points before and after the lamellae. In addition, a critical point was detected at the boundary. The location varied in different scenarios, separating two main flow zones: the critical point located at the boundary and the separation line connecting one of the lamella’s stagnation points. In the zone near the entrance, the flow reversed, and most of it returned to the inlet zone, creating a recirculation flow and increasing the average residence time. When entering the zone near the exit section, water particles could leave in a shorter time (see Figure 4).

Based on the data, two phenomena could be observed. One was the change in the location of the separation line between the two main zones, which varied with the load. The other was the presence of a swirling zone that occurred in the initial stage.

Knowing the flow pattern, the velocity derivatives were calculated in the simulation cell near the critical point, and then these were categorized based on the eigenvalues of the velocity derivative tensor for each scenario. The results are summarized in

Table 5. Signs of the eigenvalues of the velocity gradient (Jacobian) matrix near the critical points.

	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19
Re1	-	-	-	-	-	-	-	-	-	-	+	+	+	+	+	+	+	+	+
Re2	-	-	-	-	-	-	-	-	-	-	+	+	+	+	+	+	+	+	-
Im	<>0	<>0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Class	AF	AF	AN	AN	AN	AN	AN	AN	AN	AN	RN	RN	RN	RN	RN	RN	RN	RN	S

.

At the bottom part of the lamellae, the stagnation points were attachment points (attracting nodes), and the top part detachment points (repelling nodes) were situated in all cases. Near the inlet (points 1 and 2), there were foci, with a swirling flow. Nr. 19 was a saddle point, which formed a separation line with different attachment points depending on the hydraulic load. Under the DWDL condition, a separation line developed between points 5 and 19, whereas under DWDA, it appeared between points 7 and 19. Similarly, under DWDP and WW, the separation line developed between points 9 and 19.

3.3. Compartmental Model and Its Application

Based on the coherent zones, the critical points delineated the compartment architecture, which was created based on the load of DWDL; the other scenarios are also marked with dotted lines in Figure 5. Nevertheless, since the flow pattern had two main compartments, this architecture could be simplified by introducing a recirculation stream (r). In different scenarios, the simplified architecture remained the same, only the material flows and the zone volumes were subject to change. The larger the extent of the first compartment and the recirculation ratio, the higher the settling surface area; the hydraulic load projected onto the entire space also varied between the two zones. Due to the longer residence time, the HLR decreased, and accordingly, a better removal efficiency could be expected.

Zone 1 could be delineated where a portion of the mass flew back, and the degree of this phenomenon is shown by the recirculation rate marked with ‘r’. The recirculation rate indicates the percentage of the liquid that remained in the swirling zone after entry. It is possible that a portion of the fluid flow still turned back after a recirculation loop. The term ‘n’ denotes the average number of turns made in the recirculation flow. The values of ‘r’ and ‘m’ were obtained from CFD simulations; for example, in the DWDL scenario, ‘r’ = 0.29 means that 29% of the incoming flow recirculated and made at least one loop; ‘m’ = 1.52 means that the particles in the fluid flow made an average of 1.45 loops. This value came from the fact that the particle made at least 1 turn for 29% of the times, at least 2 turns for 9% of the times, and at least 3 turns for 4% of the times. In other words, precisely 20 particles made 1 turn, 5 made 2 turns, and 4 made 3 turns; therefore, the average turn was their weighted average, which was 1.45. The A_I/A ratio showed how much settling surface was present in Zone 1 compared to the entire basin. Similarly, the values of ‘r’, ‘m’, and A_I/A could be determined for other hydraulic loads. Based on this, the HLR value in Zone 1 could be calculated as follows:

$$HLR-Zone 1=\frac{\left(1+r\right)Q}{m{A}_I}=\frac{\left(1+r\right)}{m\frac{A_I}{A}}·\frac{Q}{A}$$

(17)

where Q/A is the design HLR.

Similarly, the HLR of for Zone 2 could be calculated:

$$HLR-Zone2=\frac{1}{\left(1-\frac{A_I}{A}\right)}·\frac{Q}{A}$$

(18)

The results are summarized in

Table 6. Hydraulic loading rates and their parameters in the different scenarios.

	Design HLR [m/h]	HLR—Zone 1 [m/h]	HLR—Zone 2 [m/h]	r [-]	m [-]	AI/A [-]
DWDL	0.5	1.33	0.75	0.29	1.45	0.33
DWDA	1	1.66	2.25	0.38	1.5	0.56
DWDP	2	2.70	7.14	0.43	1.47	0.72
WW	4	5.18	18.00	0.45	1.44	0.78

.

As Table 6 shows, both Zone 1 and Zone 2 received hydraulic surface loads higher than the design HLR, meaning that a higher flow rate per unit settling surface was received due to internal recirculation. In Zone 1, the largest deviation from the design value occurred at low flow rates, while in Zone 2, it occurred under the largest load. The reason for this could be that in the latter case, almost the entire incoming Q flow rate exited in the zone above and below the last lamella.

3.4. Dynamic Behavior of the Compartments

The previous results assumed steady-state conditions for the selected hydraulic loads. In order to examine the growth of the recirculation zone with increasing flow rates, dynamic simulations were required. Starting from the DWDL state, the flow rate was gradually increased in increments of 0.05 m/h from the design HLR. When convergence was reached, the proportion of the total tank volume occupied by Zone 1 was analyzed. The calculation was continued until the WW state was reached. The results are shown in Figure 6. The results indicated that the proportion of the recirculation zone gradually increased in the range from 0.5 to 1.85 times with respect to the design HLR, depending on how many lamellae were in the given zone. This might be explained by the Coanda effect, which means that the flow followed the surface of the lamellae. Above the value of 1.85 times the design HLR, the stepwise increase ceased, as the recirculation zone already occupied more than 70% of the volume.

4. Discussion

Generally, compartments are constructed based on homogeneous velocity zones, and the primary criterion for identifying the boundary between compartments is an abrupt change in the gradient of the velocity field. Nonetheless, in many cases, the definition of compartments involves the study of the velocity field in conjunction with another parameter [32]. In this study, stagnation point and coherent flow structure detection were analyzed with the critical point theory, but other techniques could be used, as reported in Alvarado et al.’s [33] study, where the zone determination was based on visual inspection, and the compartments were separated applying RTD analyses.

According to the results, the HLR design did not take into account the hydraulic phenomena within the tank, particularly the topological characteristics of fluid flow, the separation, and recirculation zones. Recirculation inside the basin could greatly influence the surface and volumetric loads, as well as the residence time. Furthermore, the appearance of recirculation could also lead to a reduction in the dead zones [34]. Recirculation is particularly useful in biological systems to maintain biomass within a system. Examples of this can be found in aquaculture [35] as well as in the field of wastewater treatment [36].

In the past, changes in flow topology were observed based on the velocity field and the related Reynolds number [37]. However, the engineering practice often works with idealized and steady-state conditions, disregarding these changes [38]. In this study, the authors aimed to draw attention to the consideration of these effects, but it was not the scope of this study to link the actual HLR to solid separation efficiency. The model did not account for the density currents that may occur due to temperature differences between the incoming flow and the temperature inside the basin [14]. Furthermore, the model considered only the settling process with discrete particles and did not include the possible flocculation. Additionally, the thickening process was also not accounted for in the model [39].

In terms of their practicality, the effectiveness of biokinetic simulation tools for water and wastewater treatment could be enhanced by integrating hydraulic phenomena, particularly in systems that contain substantial recirculation loops.

5. Conclusions

A numerical flow simulation of a lamella settler based on the finite volume method was conducted, and a topological analysis of the resulting flow pattern was performed. Compartmentalization was carried out based on the structure of the flow pattern, and hydraulic surface loading was examined in each zone. The novel technique of compartmentalization was derived from the topology of the local flow pattern. By varying the properties of the points defining the topology as a function of the hydraulic load, the extent of the compartments and the flow between them adaptively changed. The results of the study can be summarized as follows:

• Based on a numerical tracer experiment, the tracer was already detectable at around 40% of the average residence time at the outlet, while more than 2.5 times the average residence time was required for 10% of the tracer to exit. The latter may indicate the presence of dead zones, while the former may indicate hydraulic short circuits.
• The discrete phase model approximation provided acceptable results compared to the Eulerian multiphase model and the measurement results.
• The consequence of the topological analysis was that two main flow zones formed in the lamella settler, with significant internal recirculation in the first zone. The boundary between the two zones varied with different loads, and the recirculation zone became more extensive with increasing loads.
• Due to the recirculation eddies that formed in the reactor tanks, the actual hydraulic load could be much higher (from 1.5 to 4.5 times) than that in individual zones; therefore the design HLR, which referred to the entire tank, did not take into account the local flow pattern.
• Depending on the load, the extent of the recirculation zone changed in a stepwise manner until the recirculation zone became so extensive that it included all the lamellae except the last one (from approximately 35 to 75% of the total volume).