Micro-macro transition for numerical simulation of submerged membrane bioreactor

Membrane bioreactors – submerged (SMBR), sidestream ones, or combination of both, have proven their efficacy in a wide application area, including wastewater treatment and integrated processes for bioproduction and separation. Membrane bioreactors – submerged (SMBR), sidestream ones, or combination of both, have proven their efficacy in a wide application area, including wastewater treatment and integrated processes for bioproduction and separation. The integration of bioreactor with membrane separation in submerged configuration is advantageous in view of gentle trans-membrane pressure operation, imposing a controlled and constant permeate flux (1, 2), fouling control through scouring by aeration and/or solid particles and of reduced energy consumption associated with the recirculation pump in the side-stream configuration (3, 4). In recent years, a sustainable interest in CFD simulations of SMBR is observed, as they are expected to be helpful in problems related to enhanced membrane separation and effective in terms of required computational resources. CFD modeling has extensively been used for investigation of the optimal membrane module position, fouling conditions and aeration mode, including large-scale bioreactors (5, 6, 7, 8, 9, 10, 11, 12). An integrated hydrodynamic CFD model (non-Newtonian rheology, porous media approach for the membrane unit) was used to study the best membrane unit location in a full-scale bioreactor (13). The potential of 3D simulations in optimizing the inner geometry (baffle angle), gas velocity and fouling resistance was revealed in view of shear stress fouling control (8). Submerged hollow fiber unit was modeled as a porous media zone and, incorporated with a 3D two-phase CFD model, was applied to SMBR (14) for wastewater treatment. The 3D two-phase Euler-Euler method with RNG k-ε model for turbulence closure was applied for cake layer fouling control in a flat sheet SMBR (15). CFD fouling investigation was reported in a pilot anaerobic reactor with vertically submerged flat sheet membranes, the sludge mixture being treated as homogeneous non-Newtonian fluid (7). The model allows cake accumulation calculation and 2D visualization of the fouling, based on the shear intensity profiles expressed on the membrane surface.

Published studies so far have shown the potential of CFD modeling realized mainly in predicting overall hydrodynamic conditions in a given MBR geometry, and less in the study of the transfer to and through the membrane.

In a number of cases the membrane is treated as an impermeable wall, the assumption being justified by the low permeates flux as compared to the liquid velocity close to the membrane (10, 11). A quantitative framework for this condition is found in the practice of wastewater treatment by SMBR (10). Covering a wide range of ultrafiltration membranes applications, the resulting ratio of suction velocity to bulk liquid velocity is usually less than 0.01%. The assumption that the permeate flux has a negligible effect on the hydrodynamic picture has not always received an experimental confirmation, as shown in (16). The use of boundary conditions for a permeable wall accounting for the fluid flux normal to the membrane may improve the simulations in cases with higher permeation rates (17) as well as regarding the effect of the resulting wall concentration on the concentration polarization. Considering the boundary layer thickness, CFD investigations are performed both with constant or developing boundary layer (18). The latter includes calculation of the flow-field evolution in time and space under simultaneous development of the fouling layer, thus leading to an improved modeling and dynamic simulation of the membrane module operation.

The effect of the membrane on the overall reactor hydrodynamics becomes more important as the scale of the integrated unit increases. The implementation of a membrane source term in the overall balance equation, applied to the computational cells adjacent to the membrane boundary, was a successful modeling approach, supported by experimental validation in the case of fluidized bed membrane reactors for selective extraction of hydrogen via membranes (19). Treating the membrane module as a single porous media by accounting for the resistance of the latter in an additional momentum source term, was used to study the effect of membrane unit location and to optimize reactor design in a full-scale MBR (13). The study points out the differences in the obtained velocity and shear stress field as compared to the results with the impermeable wall approach. Lower, but much more uniform velocity field both in the membrane unit region and the overall reactor was observed with the integrated model.

Shear stresses relation to flux and mass transfer characteristics at the membrane surface has been reported both for cross-flow filtration cells (20, 21) and SMBR. Examples for numerical simulation of the shear field at the surface of the submerged hollow fiber or tubular membranes, coupled with an experimental investigation by the electrochemical shear probe method are given in (11, 22). Optimization of the shear stress distribution and critical particle diameter deposition on the membrane surface was searched in a CFD study of a flat-sheet sMBR in view of fouling control (15). In (23 a full-scale flat sheet SMBRs was simulated to point out the advantages of sparger bottom placement on the cross-flow velocities and shear stresses near the membrane and the homogeneity of the bubble distribution. CFD modeling of MBRs allowed calculations of the appropriate shear on the membrane surface in order to control sludge filterability and fouling (24, 25), related to sludge particle distribution (aggregation and breakage). The coupling of fouling models to mass transfer is still considered to be a challenge in view of its predictive reliability (17). The existing examples are focused on membrane elements; integrated mathematical models for MBRs have to be developed (26). In the order of increasing scale of the studies (e.g. membrane surface element, membrane module, a reactor with submerged membrane module) the number of reported investigations decreases. The analyses presented are based mostly on generalized analytical expressions, while numerical simulations are still rare. This is partly due to the large difference in the size of the membrane elements with respect to the entire reactor. One successful approach to deal with this difficulty is to model the membrane as a single porous element. This article presents a way to determine the characteristics of a porous cell based on a detailed analysis of the micro-level membrane structure.

The objective of this work is a numerical evaluation of the permeability of the membrane by simulation at a micro level, e.g. modeling of a small fragment as a representative element of the membrane structure. Then the results obtained can be transferred to the entire membrane module at the macro level.

The proposed micro-macro approach is considered to be the direct method as follows:

For our numerical experiment, we used the membrane described in (27). Schematic view and SEM micrographs of the front and backside of a polyimide microsieve with a pore size of 4 μm, the pitch of 10 μm and a grid of support bars are given on Figure 1. The membrane thickness is 1 μm.

Figure 1

The complex geometry of the membrane in the presence of a large number of pores and supporting structure leads to the need for a huge resource for direct numerical simulation. Attempting to simulate flow across the entire membrane has proved impossible with conventional computing resources (6-core Xeon E-5 with 128GB RAM for example).

This experiment gives us the idea to try to model the membrane not with its real geometry, but as a homogeneous porous medium with appropriate permeability.

The homogenization was performed by isolating a single cell at a micro-level. Here the single cell is this part of the material, by multiplication of which the geometry of the entire study area could be described. A part of the fluid in front of and behind the membrane must be attached to the cell for a description of the fluid flow through the membrane during the numerical simulation.

As it is suggested that single cells may influence each other, a representative area that contains one, two or more single cells must be identified for the purpose of the study.

Since the representative area contains few geometric elements, it can be modeled with its real geometry and successfully accomplish numerical simulation to determine the fluid flow rate across the membrane. After that, this flow rate can be used in the Darcy equation (1) from which the permeability can be determined. The calculated value of permeability can be inserted as a characteristic of the porous medium at the required place in the numerical simulation software.

Where Q is total discharge (m³/s), k- intrinsic permeability of the medium (m²), A - cross-sectional area to flow (m²), (p_b − p_a) - total pressure drop (pascals), μ - dynamic viscosity (Pa·s), L- length over which the pressure drop is taking place (m).

Methods

The simulation was carried out in a bubble column shown in Figure 2, which consists vertical cylindrical column with height 1.3 m and a diameter 1 m The sparger, located at the bottom of the column, has a toroidal shape with a ring radius 0.145 m and a tube diameter 0.120 m. Due to the axial symmetry of the object, only a 30° sector of the column is simulated, which reduced the calculation time. The main calculation area (the area over the sparger) is meshed by regular hexagonal dominated mesh.

Figure 2

Thе procedure of homogenization is demonstrated by numerical simulation of a 1 μm membrane described in (27). The space enclosed between four intersecting beams of the pad is accepted as a single cell. The geometric model of the cell, created with a CAD modeler, is presented in Figure 3.

Figure 3

It is of great importance to determine the correct boundary conditions in numerical simulation. In the case at hand, the boundary conditions reflect the interface of the insulated cell with the other cells of the membrane. Thus, they are not the product of natural physical processes but are imaginary connections that depend on the local equilibrium of a single cell. In order to reduce the effect of boundary conditions on the numerical results, a representative computational area containing more than one cell (Fig. 4) may be included in the calculation procedure. The dimensions of this representative area depend on the specific geometric structure of the cell. Prior to the basic calculations, these dimensions are to be determined in the way that the boundaries to be sufficiently distant and do not affect the results.

Figure 4

The CFD code ANSYS CFX is used to simulate fluid flows through the membrane. For this purpose, a dynamic mathematical model is developed by combining the Navier-Stokes equation and Darcy’s equation for flow through porous media which account for the membrane configuration. For the numerical solution of the equations, the second-order backward Euler transient scheme is applied. The numerical simulation is done using k-ε model with scalable wall function for turbulence, high-resolution advection scheme and second-order backward Euler transient scheme. In general, the procedure follows that, described in (28)

The software computational algorithm describes the hydrodynamics of membrane processes based on the momentum source (S_M) through the porous region using the following equations (28):

where and are the streamwise and transverse permeabilities, and - the streamwise and transverse loss coefficients, μ - dynamic viscosity, ρ – density, U – vector of velocity, U – velocity magnitude, indices x, y and z denote Cartesian components.

It is difficult to make a theoretical connection between the permeability coefficient used in Darcy’s equation (1) and momentum equations (2). Therefore, this parameter has to be initially identified. For this purpose, the water flow at a pressure drop of 5 KPa through a completely porous area 0.1 x 0.1 x1 m has been analyzed at different permeability coefficients. With these characteristics, the flow rate has been calculated once with the Darcy equation (1) and the second time using the software. The results are presented in Table 1.

Using the scheme mentioned above and the presented membrane geometry, calculations were made for the membrane permeability of different fluids at a constant pressure drop of 1 bar and variable viscosity depending on the fluid used. The results are shown in Figure 5.

Figure 5

With a constant coefficient of permeability, the flow through the membrane depends on the dynamic viscosity of the fluid. To investigate this phenomenon, a numerical simulation of imaginary fluid flow was performed through the membrane under consideration, modeled as a porous solid with a coefficient of permeability of 10^-12 m². The results for the flow rate at viscosity in the range 4x10^-4 to 4x10^-3 in Pa.s are shown in Figure 6.

Figure 6

Figure 7 shows some general results regarding the hydrodynamics of the reactor volume.

Figure 7

The almost identical results in the two columns on the right side of Table 1 indicate that the permeability coefficient of the Darcy equation is analogous to that coefficient in equations (2). The difference in values in the last row of the table can be explained by the formation of large flow velocity in the pores with greater permeability. As attempts to adjust the results with the refining of the mesh fail, this can be achieved by changing the turbulence numerics to High-Resolution advection or switch to transient analysis.

The results of the numerical simulation performed with a representative region of one, two and four cells, shown in Figure 5, show that only one cell can be used for the particular structure of the membrane under study. For all four fluids studied, the permeability values did not change significantly with the change in cell number. However, for other types of membrane structures, it is necessary to carry out such an analysis again, as their results may be different from those.

One way to separate the materials into permeate and retentate in membrane modeling with the help of the porous membrane is to define the flow as inhomogeneous and to set its own permeability coefficient for each material. However, not every software allows such a procedure. If the membrane structure shows a result similar to this in Figure 6, an additional tool may be the dynamic viscosity. Each material can be assigned an additional effective viscosity located on the steep part of the curve for permeate and on the slopes part for the retentate. The results for water volume fraction, gas holdup, water, and air velocity, obtained by numerical simulation (Fig. 7), are expected and are similar to those presented in (28).

In the numerical simulation, membrane structures can be modeled as a porous structure with an appropriate coefficient of permeability. Whether a material is permeated or retentate is determined by its viscosity and the membrane permeability

The permeability of the membrane as a porous media can be determined on the basis of a representative area (a number of membrane cells, consisting membrane, support structure and added volume of the fluid in front and behind the membrane)

The modeling of membrane structures as a porous media on the basis of a representative area makes it possible to perform the calculations on a conventional computer.