Computational Fluid Dynamics (CFD) Modelling of Porous, Ultrafiltration Membranes

This paper presents a finite element, three dimensional numerical model of flow in porous ceramic ultrafiltration membrane system together with experimental validation. The modelling is based on the computational fluid dynamics (CFD) technique. Major difficulty that arises during CFD modelling is appropriate and precise description of the porous media in terms of Navier-Stokes equation. Based on pressure-flow experimental measurements we present the approach for calculating the essential components of porous media flow resistance which are necessary for proper description of membrane process. The detailed description of a membrane was applied which accounts for support and membrane layer respectively. Moreover own procedure applied is presented, that is written using C language, for calculation of flow parameters. The model presented is in very good accordance with experimental results.


INTRODUCTION
The membrane filtration processes have played important role in the industrial separation processes. Although these processes are known for about twenty years, there are still many studies that focus on the testing of new membranes' materials and determining of conditions for optimal selectivity i.e. the optimum transmembrane pressure or permeate flux to minimize fouling [1]. Computational fluid dynamics technique may provide a lot of interesting information for the development of the membrane processes.
Ultrafiltration is a pressure-driven purification process in which water and low molecular weight substances permeate a membrane while particles, colloids, proteins, bacteria, macromolecules and other organic molecules larger than 0.01 micron size are filtered. The primary removal mechanism is size exclusion. In the ultrafiltration the porous membrane are used. Pore sizes for ultrafiltration membranes range between 0.001 and 0.1 micron. However, it is more customary to categorize membranes by molecular-weight cut off. The most popular membranes' materials are polymers (polysulfone and cellulose acetate are the most common). Due to the variety of materials and sufficient selectivity effects during the transport the polymer membranes were unrivalled until the inorganic membranes, including ceramic appeared. Since then their popularity grows, the ranks of tests are conducted to demonstrate the *Address correspondence to this author at the Institute of Technology and Chemical Engineering, Pozna University of Technology, pl. M. Sk odowskiej-Curie 2, 60-965 Pozna , Poland; Tel: +48616653771; Fax: +48616653649; E-mail: maciej.staszak@put.poznan.pl effectiveness of their use in various fields. Ceramic membranes have lot of advantages: (i) autoclavable, (ii) able to be sterilized by superheated water, steam or oxidizing agents, (iii) high temperatures resistant, (iv) acids and bases resistant, (v) solvents resistant, (vi) excellent mechanical resistant, (vii) long working life, (viii) environmental friendliness. Despite many advantages the ceramic membranes have also disadvantages -large weight and high production cost of ceramic components. However, they are compensated by the long lifespan and high chemical and thermal resistance of the membrane, especially important in the case of application uses.
There are some interesting modifications of the ultrafiltration process such as polymer enhanced ultrafiltration PEUF [2] or micellar enhanced ultrafiltration MEUF [3]. MEUF as a surfactant-based separation process is an effective technique to remove almost all the toxic metal ions and/or soluble organic solutes from aqueous solutions. In this process, the surfactant, at a concentration higher than its critical micelle concentration (CMC), is added to the aqueous solution containing metal ions and/or organic solutes. The metal ions are adsorbed on the surface of the oppositely charged micelles by electrostatic attraction. The organic solutes are solubilised in the micelles interior by ion-dipole interaction. Then the micellar solution passes through an ultrafiltration membrane with a small enough pore size to reject the micelles containing the contaminants. Micelles, the adsorbed metal ions and the solubilised organic solutes are rejected.
The permeate contains very low concentrations of un-adsorbed metal ions or unsolubilised organic solutes and surfactant monomers. Addition of surfactants, of course, changes the nature of transport through the membrane, so the modeling of fluxes for surfactant s' solution is very important due to application possibility. In the membrane filtration processes both free and porous flow occurs. The most popular equation, which describes the fluid flow is Navier-Stokes equations [4]: Where -density, v -velocity vector, p -pressure, g -gravitational acceleration vector, S -body forces (source term) and T -viscous stress tensor.
Although known for much over hundred years the Navier-Stokes equation still constitutes many difficulties when solving for the fluid velocity field v. The first two terms define the inertial term of a fluid control volume which is related to the divergence of stress and other fluid forces F (such as gravity, surface tension etc.) acting on that volume. The stress tensor (representing stresses acting on a fluid volume) is typically presented as sum of stresses acting uniformly (axiator -negative pressure gradient -p) and stresses acting in nonuniform manner (deviator -deviatoric stress tensor ·T) on a fluid volume. The most common method used to solve this vector differential equation is the finite volume method FVM known as computational fluid dynamics technique (CFD). The space of simulation is built in the form of control volumes mesh of which every element can be a geometrical entity of any type. Most often it is a triangle or tetragon and tetrahedron or hexahedron in two or three dimensional space accordingly. In general the finite element method lets to approximate a field variable (velocities, concentrations etc.) using some arbitrary basis function. This lets to represent a field variable value as a function of space location over a finite element. For fluid flow computations this function is kept possibly most simply due to a computational cost. It is typical to choose f basis = 1 as a basis function [4], although there exist higher order implementations [5][6][7]. Such choice of simplest possible basis function reduces the representation of a field variable to constant value over single control volume, typically located at its centre of gravity. On the other hand we have the spatial and temporal discretization of Navier-Stokes vector equation. Several schemes are investigated and found suitable for different purposes. The derivatives are represented as finite differences converting the partial differential equations set into ordinary differential ones in the case of dynamic simulation. If the higher order discretization is used then more accurate solution should be obtained but it is occupied by higher numerical cost and sometimes causes unstable behavior of a calculation process. Depending on the problem the first or second order spatial discretization is chosen most often. Some special techniques exist also: for example MUSCL (Monotone Upstreamcentered Schemes for Conservation Laws) that provides more sharp reproduction of specific phenomena e.g. sound shock waves.
Despite high computational demand the computational fluid dynamics is one of the best available techniques used for modeling flow in the membrane. A review of a literature revealed that most of the numerical works dealing with the modelling of a flow across a membrane, were not experimentally validated [8][9][10][11][12][13][14], or used polymeric membranes in experimental work [15][16][17][18][19][20]. To the author's knowledge no information has been reported in open literature on computational fluid dynamics (CFD) modelling of ceramic membranes. Moreover only a two-dimensional mathematical model of membrane permeation and separation was generally described.
A typical parameters of a membranes given by their manufacturers are cut-off and material from which the membranes are made. The aim of this work is to present a method, which based on basic membrane pressure-flow measurements, is suitable for accurate estimation of the resistance parameters which are needed for Navier-Stokes equations solving. For that purpose a three dimensional space CFD modelling of permeate flux for ceramic membranes and their support alone was solved and compared with the corresponding ultrafiltration process. We show also the ability to perform very fast calculation using this technique for different transmembrane pressures. The model presented is based on the finite volume method with circular symmetry condition. To validate the calculations the model results presented are compared to experimental data.

EXPERIMENTAL PART
The laboratory-scale ultrafiltration system (SPIRALB from TAMI Industries) shown in Figure 1 was used. In each experiment, the volume of the feed solution (water or aqueous solution of sodium dodecyl sulphate (SDS), at the concentration of 0.5, 1 or 5 CMC) for ultrafiltration was 1.5 L. The solution was ultrafiltered after it was adequately mixed. The experiments were performed with recycling of the retentate into the feed vessel and collecting the permeate solution into permeate test-tube. The process was run from 5 to about 10 min, depending on the transmembrane pressure. The transmembrane pressure was changed from 0.05 to 0.4 MPa. The ultrafiltration experiments were conducted at room temperature (around 25 o C).
Ceramic disc membrane from TAMI Industries with cut off 15 kDa, moreover a support alone was used in all experiments.
Membranes used in the module are in the form of ceramic discs with a diameter of 90 mm, thickness 2.85 mm and the effective surface 0.00635 m 2 . Carrier layer (0.15 mm) is made of a mixture of oxides of aluminum, zirconium and titanium, the active layer of zirconium oxide.
Anionic surfactant, sodium dodecyl sulphate (SDS) (CH 3 (CH 2 ) 11 OSO 3 Na) with a purity of 99% was obtained from Aldrich. The surfactant parameters: viscosity ( ), density ( ) and critical micellar concentration (CMC) are shown in Table 1 [21]. Measurements of dynamic viscosity of solution of surfactant were performed using a Haake's RheoStress viscometer. Density measurements were carried out using solution of surfactant on a densimeter DMA 5000, the Austrian company Anton Paar. The equilibrium surface tension was determined by the du Noüy ring's method using a K12 Tensiometer (Krüss). On the basis of these results the values of CMC of surfactants were estimated according to graphical methods as interception points of the straight lines just before and after CMC [22] (the measurements not presented in this work).

Calculation
The calculations were performed using commercial Ansys/Fluent software. The use of such a tool has advantage that it has been thoroughly tested and it allows for extending its capabilities by the use of own models. The disadvantage might be that the core of a computation cannot be modified in any manner but the need for this is a rare case.
During a flow in a porous media the significant momentum losses arise that must be accounted for. This is done by incorporating two resistance components, namely inertial and viscous ones. The inertial resistance (first term in eq. (2)), that arises due to the fluid-porous solid interaction is described by so called inertial loss term and viscous resistance (second term in eq. (2)) by the Darcy's equation. Assuming isotropy of a membrane the total momentum loss is then given by: The above equation is presented in scalar form, the subscript i denotes Cartesian direction (x, y or z), is the Darcy permeability and R is the inertial resistance factor. This equation is treated as a source term S in Navier-Stokes equation (1).
Inertial resistance coefficient R is calculated based on membrane permeate flux J measurements performed by changing a pressure drop P.
The determination of coefficient is done by relation: where L is the membrane thickness.
Due to isotropy of the membrane both resistance attain the same values for each Cartesian direction but formally they are the vector quantity. To relate both inertial R and viscous resistance parameters for membrane and support with transmembrane pressure an arbitrary, nonlinear fitting function was chosen and minimised based on obtained experimental data: The presented function f( P) is the same for both resistance parameters and R because their dependence on the pressure drop shows similar, exponential character. Based on typical least squares  The typical relations are presented on the Figure 2 where the fitted functions f( P) are drawn according to the values of a, b and c coefficients (equation (5)). The resistance of the membrane layer is almost two orders higher than support as it was expected. The values of regressed parameters are given in the Tables 2 and 3.
The three dimensional space of simulation was represented in the form of control cells mesh. The force, mass and energy balances equations were discretized over those elements constituting differential equations set which was solved by appropriate, iterative algorithm. A fundamental aspect of solving Navier-Stokes equation is the proper description of space boundaries. For the case presented the boundary conditions applied to the model were: pressure inlet establishing a transmembrane pressure and pressure outlet located opposite side of the membrane. The last boundary condition applied needs explanation: due to the expected flat velocity profile across the membrane it was decided to examine only a small, cylindrical part of the membrane by simulation. Such approach drastically reduces time of computations without losing generality of description. As a boundary condition placed at the cutting plane, the so called symmetry boundary condition was chosen. That particular kind of boundary condition represents boundaries of a part of a membrane, assuming zero normal velocity and zero normal gradients of all variables at a symmetry plane.
The diameter of the circular part of a membrane was chosen to be 1mm because it is big enough to characterize the flow. The following parameters (size of mesh) was applied: total number of tetrahedral cells = 54942, divided into two zones: 26647 tetrahedral cells in membrane layer, 28295 tetrahedral cells in support layer (Figure 3). Because the Ansys/Fluent 14.0 software allows to extend the main computational core by attaching own models and algorithms at many different points of computation process, it opens many research possibilities to enrich the CFD modelling technique. The resistance of a matrix and membrane layer was defined by formulating it in terms of User Defined Function (UDF). That UDF was created using C programming language and compiled in the Ansys/Fluent computing environment. The formulation was created by the use of DEFINE_PROFILE subroutine and attached to volumetric cell threads representing matrix and membrane, respectively. Additionally the integrating procedure was created to compute the volumetric mean integral of permeate flux J according to equation (6): The integration was performed at the end of every succesful iteration and time step computation (DEFINE_EXECUTE_AT_END). The used routines of UDFs are presented in Appendix. The model solver chosen is unsteady formulation and supply pressure was changed linerarly from 0.05 to 0.4 MPa during 1.0 seconds of simulated time. It is assumed that achieving a momentum steady state is very fast process especially in the case of noncompressible flow. The time step was arbitrarily chosen as 0.01 s. The convergence criteria for time step numerical acceptance was to achievie the value of the residuals of the solution less than 0.001. That means that obtaining less than 0.1% of solution imbalance is accepted in the solution.
The porous membrane and support zones were characterized by inertial resistance coefficient in accordance with direction vectors. In our case we assumed isotropy of the support and membrane, but such vector definition of inertial forces allows also the anisotropy of the porous media to be accounted for.

Experimental Validation
Employing Darcy's law permits the permeate velocity to be expressed as a function of the pressure difference existing across the membrane wall. Determining the pressure along the surface of the membrane allows a straightforward calculation of the permeate velocity profile along the membrane wall. The pressure outside the membrane module is assumed atmospheric. Figure 5 shows the gauge pressure profile and velocity vectors. In the first stage by smoothly changing the transmembrane pressure range between 0.05 to 0.4 MPa the distribution of pressure across the membrane was calculated. Moreover the velocity vectors were obtained. As expected, an increase in transmembrane pressure caused an increase in transport of water across the membrane filtered.  Because there is no possibility of experimental verification of pressure distribution on the membrane, in the next figures (Figure 6) only a comparison of filtered water velocity, the permeate flux depending on the transmembrane pressure applied are presented, for water and the solutions of SDS, respectively. As it is shown the model describes well the experimental data, the calculated deviation of residual indicated small differences between experimental values and the values obtained from CFD calculations.
The study of the relative model deviations (presented in Table 4) gives the view that the obtained calculation results are accurate with the error of magnitude of a view percents. Only one case exceed 50% of relative error but that seems to be accidental measurement error, that is clearly visible as first experimental point on the Figure 6c. This leads to the conclusion that the proposed, viscous and inertial resistance components calculation approach for the porous media, and consequently the CFD technique itself, are in very good accordance with experiments.

CONCLUSIONS
The method for estimation fundamental resistance parameters for transmembrane flow is presented. To validate the appropriateness and accuracy of presented approach we solve the Navier-Stokes equation for flow in porous media. The three dimensional transient finite element model for the prediction of flux in the ultrafiltration system is presented in this study. The CFD model was setup using tetrahedral mesh representing two layers: membrane and support. The Navier-Stokes in three dimensional space were discretized on the mesh control volumes constituting ordinary differential equations set. The nonlinear relation between transmembrane pressure and resistance in the porous media was also introduced into the model. The solution obtained shows typical flow behaviour -velocity along both membrane and support zones is constant, while the total pressure drop occurs mainly at the membrane layer which is characterized by a resistance of two orders of magnitude higher than support zone.
The future development of proposed method will consist on adding a model to relate the resistance parameter with the pressure drop. This approach applies an arbitrary exponential equation in that part of computation which for some rare cases might sometimes reveal inaccuracies.  C_PROFILE(c1,t,i) = (a / (b*pressure + c)) + d;

SYMBOLS
end_c_loop(c1,t) } The function above formulates the inertial resistance term by relating it to the transmembrane pressure according to equation (3). The conditional block is used to achieve a linear transition region, from the start to the first second of the process, which results in faster convergence to the solution. This formulation is in accordance with pressure inlet boundary condition. The exact attachment for the cells is formulated as a inertial resistance term in momentum equation (1).
A function that computes the permeate flux J according to equation (6)