Modeling and simulation of membrane process

2.1.1 Standard blocking model

The model is also referred as pore blocking model. The rate of change of number of pores that are available for fluid transportation at any time of the process can be given by relation:

In the above n represents number of non-blocked pores accessible by the fluid per unit of area of the membrane surface. The variable f_b is the number of pores that can be blocked by a unit volume of particles. The C₀ is the initial particle concentration in the feed stream. The indices 0 refer to the initial values of specified variables. The analytical solution to the above is:

Thus two time-dependent component fluxes can be identified, the flux by unblocked pores and flux by blocked pores. The dynamics of the total flux can thus be given:

where u_b0 is the flux for the blocked membrane.

2.1.2 Cake formation model

For incompressible cakes the Carman-Kozeny [4] equation relates pressure drop with velocity and membrane details:

Equation commonly used to estimate specific cake resistance r_c reads:

In the above equation, μ is the viscosity of the fluid, ε is the porosity of the solid cake, D_p diameter of the particle, u_s superficial velocity, Φ_s is the sphericity of the particles in the cake. The time required for 99 % reduction of the initial flux is given by:

where c_g is the particle concentration on the surface of the cake layer, c₀ is the particle concentration in the bulk (feed) fluid, R_bm is the resistance of the blocked membrane and u₀ is the initial velocity (to the clean surface of the membrane).

2.1.3 Physicochemical model

The model with a fouling limitation applied is given in the form that reads [5]:

where R^* is the long-term fouling resistance, p, q, and r are specific constants. The variable p renders very fast initial increase in the resistance, due to the substance adsorption at the surface of the membrane which causes blocking of the pores. The q and r variables render the curvature of the time evolution of the hydraulic resistance of the membrane. R^* is given by a Langmuir or Freundlich isotherms and is a function of bulk concentration and solution pH. Equation 13 represents typical correlative approach where the specific variables must be experimentally obtained for individual membrane type.

2.1.4 Reaction model

The assumption of Langmuir–Hinshelwood kinetics for the buildup of an adsorbed layer on the membrane surface [6] and leads to the description of the dynamics of its thickness l:

The C_w is the concentration of the separated substance at the membrane surface. The kinetic constants k₁ and k₂ refer to the rate of its deposition and removal respectively. The solution to the above reads:

The values of kinetic constants can be obtained by means of nonlinear fit numerical methods versus experimental data of cake thickness formation.

By solving classical diffusion equation of Fick’s second law, other formulation can also be given:

where ΔP is the pressure across the membrane, μ is the viscosity, p is permeability of the filtrated substance, l represents thickness of the layer, C_b bulk concentration of the substance, and R_m is the membrane resistance. The kinetic constants represented by k_r and k_m refer to the deposition and removal process respectively. The solution to the above problem can only be done on the numerical basis.

2.1.5 Maxwell permittivity model

Gas separation by selective transport through mixed matrix polymeric membranes is commonly used technique for gaseous mixtures treatment. Chung et al. [7] suggested to use Maxwell permittivity model [8] originally aimed to predict the material resistance to external electric field. The effective permeability of composite membrane can be formulated in terms of component permeabilities of continuous P_c and dispersed phases P_d:

The above model is a simple approach that accounts for dilute suspension of spherical particles. The variable ϕ_d is the dispersed (polymer) phase volume fraction. Some modifications of the above approach have been proposed [9] to account for more realistic behavior. The diffusion coefficient for gas in membrane pores is then given by Knudsen diffusion coefficient that reads:

where r is the effective pore radius expressed in Å, T is the absolute temperature, M_A is the molecular mass of the gas, and d_g is the diameter of the gas molecules also expressed in Å. The value of the permeability of the combined sieve (P_d) and interface (P_i) void can be used along with the continuous polymer phase permeability (P_c) to obtain a predicted permeability P_3MM for the three-phase mixed matrix materials by applying the Maxwell equation:

The variable ϕ_d is the dispersed (polymer) phase volume fraction and ϕ_i is the volume fraction of the interface void.

2.1.6 Modeling dense membranes

The nonporous membranes, i.e. for reverse osmosis [10], modeling shares the same principia of mass conservation as the porous, they differ in the form of constitutive relations although. The key process governing the transport through the membrane material is the solution-diffusion mechanism. The performance of the process depends on the selectivity φ_M of the membrane which is defined by:

where C_w and C_p are the concentration of diluted substance in the feed and permeate, respectively. The flux of permeate J_P can be calculated by relationship:

where the thermodynamic function to compute a maximum value of the selectivity coefficient of the membrane reads:

In eq. 21 the variable j_m is the impermeability coefficient of the membrane, and ΔP_m is the operating pressure. The ω is an empirical coefficient that depends on the physical properties of the material and the operation conditions. The ΔW is the component of power of osmotic process expended for transporting ions across the membrane, T is the temperature of the fluid, and k_B is the Boltzmann constant. The flow through the membrane under the operating pressure ΔP_m, including impurities, depends on the hydraulic resistances of the membrane R_m and the deposit layer (cake) R_c:

where the d_c and d_m are the cake and membrane thicknesses, μ_c is the dynamic viscosity of the deposit layer.

The steady state and dynamics for dense membranes can also be modeled by considering that the driving force of the process is entirely of diffusive nature [11]. In this case the flux of specie i^th is related with diffusion coefficient and specific formulation of driving force which includes corresponding permeabilities p_i,perm and p_i,feed difference:

where D_i is the diffusion coefficient of specie i^th, S_i is the sorption coefficient, H_i is Henry constant, and L is the membrane thickness. The postulated dynamic description relies on variation of diffusion coefficient and permeabilities by introducing their functional forms dependent on time, D_i(t) and p_i,perm(t). Such mathematical simplification is valid only with assumption that diffusion coefficient is independent of diffusing specie concentration which is valid for dilute solutions. This statement is postulated by the authors [11] exclusively by particular formulation of Fick’s second law for time evolution of specie concentration in membrane material:

The initial condition for the above reads:

which formulates that no diffusing specie is contained inside the membrane in the initial moment of the process. The boundary conditions read:

{c(0,t)=c1−D(t)∂c(L,t)∂x=J(t)

that depicts that at the membrane surface (x=0) there is some initial concentration c₁ of the solute, and at the downstream surface of the membrane there is a flux J(t) given by the partial pressure in the permeate channel during the process.

2.1.7 Models dependent on membrane configuration

2.1.8 Summary

The presented section outlines the broad spectrum of model approaches that are based on traditional engineers and designers techniques of membrane processes description. The modeling and the designing procedure are the tasks that employ the same tools but take different view of the process considered. The correlations presented can be used to make a simulation if some initial data like feed source, feed quality, feed/product flow, and required product quality are given. Then having selected the flow configurations and membrane elements type and setup, one can determine whether the process works as expected. The design procedure on the other hand implicates a specific results to be obtained, namely the material and geometry estimation of the membrane itself. Both views, the modeling and the design, are essentially the same regarding the theoretical and mathematical description of the process but target different needs. Most of the presented formulations are based on correlated experimental results applied to equation of mass or momentum balances. In general, the fundamental disadvantage of the models based on experimental measurements or any specific material constants obtained under particular conditions is their higher, compared to other approaches, tendency to decrease the accuracy and/or reliability when performing calculations for systems that differ significantly from systems for which they were originally established. Many correlations are designed for some specific range of temperatures, pressures, and so forth and beyond these ranges they might not be reliable source of knowledge about the system. The exception from such typical correlative issues might be to rely only on fundamental physicochemical data which do not depend on the system conditions (e.g. critical pressures, temperatures, or volumes of species involved). The mathematical approaches and models presented in next parts of this chapter address this problem because their fundamental approach is not to rely on actual, narrow range specific parameters.

2.2.1 Quantum mechanics models in polymer membranes

The underlying fundament of all quantum methods is Schrödinger equation which defines molecules in terms of interactions among nuclei and electrons and depicts geometry of the molecules by structure satisfying minimum energy requirement. The time dependent Schrödinger equation reads:

In the above the Ψ represents wave function representing quantum state of the molecule, atomic, or subatomic system. A single wave function covers at once all the particles in the system analyzed. There are many different wave functions for single system though, each of them representing appropriate description, in a chosen coordinate system. Every of the distinct representatives express all observable elements of information about the state (e.g. position, translational momentum, orbital angular momentum, spin, total angular momentum, energy). Ĥ is the Hamiltonian operator that acts on a function Ψ. The Hamiltonian corresponds to the kinetic and potential energies of a system and is the summation of all energies for every particle in the given spatial configuration at an instant time. The Schrödinger equation cannot be solved directly and approximations must be used in order to obtain results. The most commonly approximation used in order of increasing simplifications applied are as follows. The Born–Oppenheimer approximation with the assumption that the atomic nuclei do not move and only electrons kinetic energy, and the nuclei–nuclei charge interactions are accounted for. The Hartree–Fock approximation is further simplification that assumes no interaction between electrons themselves. This leads to formulation of atomic (and molecular) orbitals which can be related to stationary state of wave function for an electron in atom or molecule. The orbitals render the structure of the atom or molecule as if all the electrons formed the molecule shape without electron–electron interactions, which obviously is a simplification. The LCAO approximation (Linear Combination of Atomic Orbitals) is used in Hartree–Fock approximation to render the linear combination of atomic orbitals to further calculate the molecular orbitals by their quantum superposition. The molecular orbitals are expressed in the sense of so-called basis sets which are the linear combinations of a set of functions. The basis set is created of a finite number of orbitals, centered at each atomic nucleus within the molecule. The Hartree–Fock approximation is applicable to molecules composed from up to hundreds of atoms. Assuming no electron–electron interactions and treating every electron as independent of each other leads to overestimating of electron repulsion energy which consequently render as too high total energy of the system calculated. To overcome this problem a correlated models were proposed in which the electron correlation accounts for coupling the electrons motions that result to decrease the repulsion energy of electron–electron interactions. Lower repulsion energy estimated leads to lowering the total energy. The density functional models propose approximate correlation by an explicit term formulated on the basis of Density Functional Theory. The main assumption of this approach is the fact that the sum of the correlation energies of uniform electron cloud can be evaluated based on its density. The density functional models are usually appropriate to estimate equilibrium, geometries of transition-states, and vibrational frequencies. The density functional models reveal some serious deficiencies with regard to energetics of reaction. Due to this problem the Configuration Interaction Models were introduced in which the “auxiliary” wavefunction terms corresponding to excited states are additionally included with regards to Hartree–Fock method. The Configuration Interaction Method approaches the exact Schrodinger equation solution if the basis set applied is more complete (covers wider range of possible wavefunction states). The importance of using more complete basis sets arises due to size consistency, where typical problem is that the solution to the energy of two infinitely separated particles is not double the energy of the single particle. Due to the high demand on computing power these methods are limited to relatively small systems, however. The simplest model offering improvement over Hartree–Fock method is Moller–Plesset theory of second order (MP2). Also higher order models MP3, MP4 were formulated but they are limited to small systems. The MP2 method is practical for molecules of moderate size, although.

The basis sets mentioned earlier are functions introduced first by Slater [17] to give a description of electron orbitals. Slater type orbital (STO-xG, where x is the number of Gaussian functions) basis sets are simplest sets of Gaussian functions representing orbitals. It is significant that the same number of Gaussian functions is used to represent core and valence orbitals. They give rough estimation for orbitals but are much less computationally demanding. Pople et al. [18] formulated more complex basis sets named split valence basis sets. This formulation represents the valence orbitals by more than one basis function due to the fact that valence orbitals take main part in molecular bonding. A lot of basis sets are discussed in the references [19, 20, 21, 22, 23].

Due to the complexity of the quantum calculation it is typically not possible to perform the calculations by homebrew computer codes, instead a commercial (Gaussian, Jaguar, etc.), academic (e.g. Gamess, Firefly) or free General Public License (CP2K, PUPIL) codes are used.

Taking into account polymer membranes the main research area where the quantum modeling calculations are mostly involved are polymeric membranes in fuel cell simulations, namely proton exchange membrane fuel cell, PEMFC (Figure 2). Fuel cell is an electrochemical device that transforms the chemical energy into the electrical energy without any assistance of additional component. Principal elements of a single PEM fuel cell are: membrane, catalyst layer, gas diffusion layer (also referred to as porous transport layer), current collectors, and flow channels. Such membranes are often denoted as polymer electrolyte membranes or proton exchange membranes. The hydrogen is the fuel for the PEMFC and the hydrogen ion, proton is the charge carrier.

Figure 2.

Fuel cell typical setup.

The anode reaction reads:

The cathode reaction reads:

The overall reaction is then given by stoichiometry:

The anode reaction reads:

The cathode reaction reads:

The overall reaction is then given by stoichiometry:

One of the most popular membrane material present on the market is Nafion, which is sulfonated tetrafluoroethylene based fluoropolymer-copolymer (RF[OCF₂CF(CF₃)₂]_nOCF₂CF₂SO₃H). Nafion is most often formed as a thin foil and used as a protons transporting membrane, namely ionomer, while being a barrier for electrons or anions. This material exhibits also high thermal and chemical resistance.

Yu et al. [24, 25] simulated chemical mechanisms of degradation of Nafion using JAGUAR 7.0 code. They incorporated Density Functional Theory to formulate the potential energy surfaces for various possible reactions accounting for intermediates which could be created when Nafion is exposed to molecular hydrogen, molecular oxygen, or protons in the presence of the platinum catalyst. The B3LYP (Becke, three-parameter, Lee-Yang-Parr) exchange-correlation functional was used with 6-311G** basis set (six Gaussian primitive function for core orbitals, and a linear combination of two basis function for each valence orbitals, the stars indicate extension to the main basis set by adding double polarization functions). Based on the results three mechanisms by which OH radical species formed on a platinum surface lead to degradation of Nafion polymer were proposed and experimentally validated by F NMR and mass spectroscopy.

Sakai et al. [26] used Gaussian software commercial tool to model the deprotonation reaction of the sulfonic group (SO₃H) in a hydrocarbon membrane, as a new proton conductor for polymer electrolyte membranes (PEM) of polymer electrolyte fuel cells (PEFC). Similarly as case previously described 6-311+G(d, p) basis set was used at the B3LYP density functional theory level. The difference in basis set is presented by the plus sign meaning that the diffuse function is used for p and d orbitals. The diffuse functions are of high significance for anions to describe the electrons that are at longer distances from the nuclei. The quantum model was used to investigate and simulate the behavior of the sulfonic group of sulfonic acid C₆H₅SO₃H as a model of a hydrocarbon membrane. As the concluding statement the following two properties were noticed as the main factors decreasing proton conductivity at low hydration levels: sulfonic acid has low acidity because deprotonation of SO₃H group and protonation of SO₃⁻ anion occur with the same probability; and the SO₃⁻ anion of sulfonic acid binds strongly with the H3O⁺ in comparison with short side chain perfluorosulfonic acid. Paddison et al. [27] presented earlier similar study using B3LYP/6-311G** level density functional theory to perform modeling of a hydrocarbon polymer membrane (CF₃CF(O(CF₂)₂SO₃H)(CF₂)_nCF(O(CF₂)₂SO₃H)CF₃, where n = 5, 7, and 9). The quantum simulation performed exposed that the number of water molecules needed to connect the sulfonic acid groups is a proportional function of the number of fluoromethylene groups in the main chain, with one, two, and three water molecules needed to connect the sulfonic acid groups in parts with n = 5, 7, and 9, respectively. It is also worth noting that the simulation was performed with explicit addition of water molecules to each of the polymeric fragments (an implicit approach would be to replace intermolecular interactions by so-called averaged field), which allowed to conclude that the minimum number of water molecules required to affect proton transfer also increases as the number of separating tetrafluoroethylene units in the main chain is increased. These quantum modeling calculations provided the results that might be used to evaluate the impact of changes of polymer chemistry on proton conduction, including the side chain length and acidic functional group.

There are also modeling studies which are based on approaches which in their background rely on quantum calculations. One of such approach is COSMO model (conductor-like screening model) which is used for estimating the electrostatic interaction between a molecule and a solvent. This model belongs to continuum solvation or averaged field group of models. The main concept is to approximate the solvent presence by a dielectric continuum, surrounding the molecules. The COSMO method can be used with Hartree–Fock methods or density functional theory calculations. The extensions for real solvents allow to estimate the chemical potentials of substances in liquids which consequently leads to other thermodynamic equilibrium properties (such as solubility, activity coefficients, vapor pressures, partition coefficients and free solvation energies) calculations. The example of the use of such an approach is the work of Shah et al. [28] They used a COSMO-SAC model (conductor-like screening model – segment activity coefficient) to estimate the structure–properties relationships between the polymer structure and membrane performance. The GGA/VWN-BP functional setting and the DNPv4.0.0 basis set was used. The generalized gradient approximation (GGA) takes into account variations in the density by including the gradient of the density in the functional. The Vosko-Wilk-Nusair [29] (VWN) is the correlation for many electron system of the spin-polarized homogeneous electron gas. The VWN correlation energy per particle of homogenous electron gas improves the theory and experiment agreement. The VWN-BP represents the Becke-Perdew version of the VWN functional. The model can successfully predict two component sorption from single component sorption data. The model functions poorly in the case of some polar solvent/polar polymer combinations, however.

2.2.2 Molecular mechanics models in polymer membranes

Molecular mechanics describes molecules and chemical bonds based on van der Waals and Coulombic interactions. This approach utilizes classical mechanics laws to model the molecular systems.

The molecular mechanics approach expresses the energetic state of a molecule in terms of a sum of contributions resulted from deviations from an ideal molecular state of lowest energy. The summation is taken over stretch contributions (bond distances), bend contributions (bond angles) and torsion contributions (torsion angles) of mechanical state of molecular system. Additional terms are the contributions due to van der Waals and Coulombic (non-bonded) interactions. Stretches and bends contributions of molecular bonds are formulated in terms of Hook’s law forms. Torsion contributions need to reflect the torsional potential periodicity which takes the form of trigonometric functionals that accounts for such behavior. The energy of torsion can be formulated in terms of Fourier series which allow to account for energy differences in cis and trans conformers or between planar and perpendicular conformers. The van der Waals and Coulombic interactions summation formulate the non-bonded interactions terms. Some explicit summation terms might also be included to account for such interactions as hydrogen bonding. The van der Waals interactions are formulated as a summation of attractive and repulsive components. The Coulombic term is formulated to take the charges interactions into account.

The advantage of molecular mechanics models over quantum approach is their greater ability to transfer the geometrical variables between molecules, and expectable parameters dependence on atomic hybridization. This makes the molecular mechanics approach to be easier when applying initial geometrical guess to a numerical solving procedure and allows to perform faster calculations. In practice the computational load increases linearly with the size of molecular system modeled (which is not the case for quantum mechanics). The parametrization of molecular models is also important advantage over quantum approach. On the other hand the simplicity of the molecular mechanics models imposes some limitation when comparing with quantum approach. The most important limitation is that molecular mechanics allow to represent only the equilibrium states (geometries and conformations) of the molecular systems. The parametrization of functionals describing force fields has not yet been done to comply with nonequilibrium molecular states. The next limitation is that the molecular mechanics cannot describe the electron distributions in molecules which is an important feature when simulating chemical reactions. Consequently, the underlying principal nature of molecular mechanics modeling is its approximation character, for which the accuracy depends on quality of the parameters used and on the classifications of relations between molecules. This is typical that for such an approach going beyond the range of parametrizations, i.e. describing some new molecular structures, is expected to be burdened with high and difficult to estimate error.

The extension to molecular mechanics approach is molecular dynamics, which besides of the molecular mechanics interactions, consist on solving, typically on numerical basis, the classical equations of motion. After the calculation of potential energy field, the forces acting on the molecules are computed. Applying the Newton’s equations of motions to the molecular system results in a molecular momenta estimation which consequently leads to movement patterns of molecules. Such a system of differential equations is solved using specified step-by-step numerical algorithm of integration. The differential equations set is typically characterized by high stiffness (relation of fast to slow dynamic components of the process simulated) and special numerical treatment is needed. Historically the Gear’s algorithm [30] was designed for such systems, which is the two step method assembled of explicit Adams and implicit backward differentiation formula. Recently Verlet algorithm [31] is used to solve this problem, which is purposely designed to integrate the equations of motion.

As an example the effects of molecular structure on the polymer membrane process is studied by means of molecular mechanics calculations by Shen et al. [32] The molecular mechanics models require to input initial structure of the modeled molecular system. For the polymers a heuristic method [33] can be used, which consist on building the structure from monomers and locating them randomly in computational space. Authors examined the rejection mechanisms for contaminants in polyamide reverse osmosis process and based on the molecular dynamics calculations concluded that the solute rejection is positively correlated with the Van der Waals size of the dehydrated solutes due to decreasing accessible intermolecular volume between polymer chains of the membrane with solute size. Devanathan et al. [34] performed an ab initio study (from first principles) of proton transport in Nafion polymer membrane. One of the most fundamental problems in proton transport is that it does not follow diffusive transport of water molecules but make use of surrounding hydrogen bonded clusters. For complex molecular systems, such as those which appear in proton exchange membranes, molecular dynamics methods that aim to explain proton transport accurately need to capture this dynamical behavior. The study was conducted for different hydration levels: dry, hydrated, and saturated fuel cell membrane. In the confined environment of the membrane, the protonic defect exists as H₅O₂⁺, H₇O₃⁺, and H₉O₄⁺ cations for the three hydration levels respectively. As the hydration levels increases, the barrier for proton transport from contact ion pair to solvent separated ion pair decreases, the water molecules form a penetrating cluster, which allows the proton to hop through the water cluster structure contacting multiple sulfonate groups of Nafion polymer. The estimated by means of molecular dynamics diffusion coefficient for proton transfer is equal to 0.9·10^–5 cm²/s, which for high hydration levels stands in accordance with experimental trends, as stated by authors.

The diffusive transport in four polymer membranes is investigated on the basis of molecular dynamics by Ling et al. [35] to give the estimation of diffusion coefficients in this process. Based on diesel components the study presented the modeling calculations that allowed to investigate the micro diffusion mechanism inside polymer membranes. Authors took into account three key factors affecting the diffusional transfer: interaction energy between molecule and polymer membrane, fractional free volume of the membrane and mobility of polymer chains. The mean square displacement formulation was used to calculate the diffusion coefficient from the results obtained:

where r_i(t) is the spatial coordinate of molecule i at time t, and r_i(0) is the initial spatial coordinate of molecule i. The diffusion coefficient is calculated from Einstein formula:

where d is the dimensionality of the Brownian movement (space dimensionality), t is the time over which the diffusion coefficient value D is averaged. The power law exponent α is in normal diffusion regime (the diffusing molecules jump between individual free spaces in the polymer membrane) equal unity. The calculated diffusion coefficients were of order 10^–9, m²/s which is in good accordance with presented experiments for materials investigated.

2.2.3 Summary

In general quantum and molecular mechanics simulations even of idealized polymer membranes can provide valuable guidance by predicting the performance of candidate systems as a function of the structural and compositional factors that influence transport properties. Results obtained from such calculations, even if semiquantitative, can be used to give directions of synthesis toward suitable materials among the increasing number of potential polymers combinations. The modeled and calculated wave function representing the state of molecule can effectively be related to chemically valuable concepts such as atoms in molecules, functional groups, bonding, the theory of Lewis pairs and the valence bond model [36].

3.1.1 Finite volume method

A practical realization of solution algorithms to the Navier–Stokes equations is a finite volume method widely known as CFD. The finite volume refers to the discretization approach used, which is most popular in fluid dynamics simulation software packages. The governing form of the finite volume equation reads:

which formulates the principia that changes in time are balanced by changes in space. The vector of conserved variables Q is integrated over a volume V, while fluxes vector F is integrated over boundary surface A surrounding fluid volume V.

The Navier–Stokes equation solution is obtained by applying appropriate discretization scheme over time, for time derivative of function u:

In this way the partial differential equation is discretized to give ordinary differential equations set. The spatial derivatives follow similar schema (although different schema are used), for the x component of u:

The above scheme represents a general guideline of discretization. In typical practice the geometry of analyzed object is constructed by spatial mesh which can be assembled of number of cells. Such volumetric cells can be chosen as tetrahedrons, wedges, prisms, hexahedrons and polyhedrons for three-dimensional cases. For modeling in two dimensions it is typical to use triangles, rectangles or quadrilateral to represent the space of the flow. It is important to mention that the geometric primitives used are of convex type.

The process of mesh creation is the preprocessing stage where a decision is made about the sizes and shapes of boundaries, eventually moving or rotating bodies, and the type and size of the mesh cells used. This is crucial to the accuracy of the solution for the cell size not to be larger than the smallest process phenomena to be tracked and simulated. From this point it is evident that to be able to solve or represent very small flow effects one should also use very small cell sizes, which in turn would increase the cells number. Because the differential equations are discretized over the cells producing corresponding number of algebraic equations, the higher number of cells results in larger equations set to be solved. Consequently as a reason of increasing demands on accuracy, decreasing size of cells produces CFD models that require more time to reach the solution. This is a problem of tradeoff between the demanded accuracy and computing load when performing the step of preprocessing. The preprocessed mesh is then a base of the core step – the processing that consists on applying the iterative algorithms to obtain the solution for steady or transient states of the flow. The final stage is to postprocess the results obtained to analyze them. It can be the graphical representation of the flows but also calculation of several process variables, e.g. calculating flowrates, pressure drops, or chemical reaction rates.

Reformulating the Navier–Stokes eq. (50) to account for transport through porous media an additional source force term S must be introduced:

which accounts for the resistance effects when modeling permeate flux. For any type of membrane, polymer, ceramic, and so forth when the sieve transport is considered the most straightforward way to calculate the S resistance term is to use Darcy equation.

and inertial losses factor:

In the above formulation the multiplication of both the vector R (inertial resistance factor which for non-isotropic media depends on direction) and vector u is a scalar product. The velocity vector magnitude is designated by |u|. The resistance term S in the Navier–Stokes equation is the summation of both component resistances. For the case when the membrane material can be assumed to be isotropic, then the geometric direction i^th component of S reads:

The variable α is the permeability and R is the inertial resistance factor. For the isotropic case the indexes are omitted. The simple and straightforward estimation of both the parameters can be done on the basis of experimental measurements of the flows and pressures across the membrane [40]. Inertial resistance term R can be calculated based on measurements of the permeate flux q subject to changing transmembrane pressure Δ p:

Having calculated the R coefficient and knowing the membrane thickness L the permeability coefficient is obtained by relation:

Many software tools aimed for CFD calculations which allow fluid transport in porous media simulation, provide specific interface to enter these component resistances, in different range of functionality, however (e.g. Ansys/Fluent allows for vector resistance specifications while CFX allows only a scalar definition, concerning version 16.0).

Pellerin et al. [41] used Darcy law to simulate a two-dimensional case of membrane flux. The fluid flow outside membrane is simulated using closure model for turbulence effects (density-weighted ensemble-averaged, which is a simpler case of Reynolds averaged stress model). The membrane itself is simulated as an isolated entity by the laminar flow approach with velocity calculated on the basis of Darcy equation. The two-dimensional space was constructed by regular mesh of squares which formulated simple discretization problem of finite differences approach. The results were compared to analytical models of Berman [42] and Gupta [43]. The model proposed and both literature models show similar results which support it accuracy.

The CFD simulation of porous membrane of fuel cell is presented by Wang et al. [44] For the first time, the polarization curve and flowing components concentration distributions encompassing both single- and two-phase regimes of the air cathode are presented. The liquid and vapor transfer is governed by capillary transport mechanism and molecular diffusion, respectively, due to negligible small air velocity within the porous electrode. Authors applied following simplifying assumptions to the CFD model presented: the gas diffusion is isotropic and homogeneous, and characterized by an effective porosity and permeability of the porous media; the catalyst layer is treated as an infinitely thin surface for oxygen reduction at the membrane interface forming water; the cell is at isothermal conditions; the gas phase is an ideal mixture. The model construction allows for algorithmic coupling to predict the transition between the single and two phase regimes of flow. The two-phase model is appropriate to model the air cathode of both hydrogen and direct-methanol fuel cells.

The three-dimensional model of polymer membrane fuel cell was presented by Siversten et al. [45] The model proposed was non-isothermal model for proton exchange membrane fuel cells, which was implemented into a CFD code. The model algorithm presented allows for parallelization of the computing architecture. The convective and diffusive transport mechanism was accounted for which allows to present the concentration distribution of components. Appropriate energy balance is used to render the distributed heat generation associated with the electrochemical reaction in the cathode and anode. Authors implemented the modeling procedure in commercial tool Fluent 6.1 by the use of so-called user-defined functions, which allows to program own algorithms using C programming language. The modeling algorithm included flow channels and membrane itself. For the membrane the electro-neutrality was assumed, which is a valid assumption for fully humidified polymer material. The heat transport through the membrane for N-component mixture is calculated by the energy balance equation that reads:

where E is the total energy, k_eff is the effective conductivity, J_j is the diffusive flux and H is the enthalpy of the j^th specie. The effective stress tensor matrix is designated by τ_eff (is the volume average of the stress matrix [46] as given by (48)) and S_h is the source term for energy per unit volume and unit time which may include the heat of reaction. A general thermodynamic approach is used to model the transport in the membrane material. The current density of protons j is given by:

where L₊₊ is the specific proton conductivity [47] (assumed at value 13.375 S/m), ϕ_m is the local potential, Δ _w is the local chemical potential of water. The water flux is given by:

The L variables are the Onsager coefficients, and can be related by L_+w=L_w+. They can be estimated by the equation:

where ξ is the number of water molecules transported with each proton (value of 3 is used). The variable F is the Faraday constant. The variable L_ww is not discussed in the paper unfortunately, the reader is referred to master thesis of Sivertsen [48]. The l_mem variable is used (permeability of the membrane, 2.35·10^–7 mol²s/m^–3kg^–1):

The number of cells of a mesh constituting the space of simulation was total of 546,000. Such a mesh resolution is typically sufficient to obtain accurate and reliable results. Global comparisons show that the model and experimental results are in good accordance. The study also demonstrates that altering the conductivity modifies the current distribution by changing the relative influence of electrode potential to activation overpotentials (measure of electrode polarization).

Wiley et al. [49] present a two-dimensional study to compare the CFD simulations results with previous literature experimental data. The developed approach was a general purpose model of concentration polarization and fluid flow in pressure driven membrane separation processes. The model treats the channels as a typical space of flow with no slip boundary conditions at the walls. The membrane wall is treated differently although, where the tangential velocity component is applied to be zero and normal component of velocity (entering the space of membrane material) u_w is defined by:

where u₀ is the velocity at channel inlet, c_w is the concentration at the wall, c₀ is the channel inlet concentration. The coefficient β is given by Brian [50] as:

and R (fractional component rejection/retention):

where J is the component flux through membrane. The variable π₀ is the osmotic pressure of the feed solution, and Δ P is transmembrane pressure.

The normal velocity u_w to the membrane surface can also be given by relation of Doshi et al. [51]:

although no further details are given referring constants B₁ and B₂. The variable h represents the half flow channel height, and Sc is the Schmidt number. The main concept is to balance the convective flux in the channel with the diffusive flux in the membrane material, taking into account that not all of the fluid pass through the membrane. This last is taken into account by R rejection coefficient. Finally the balance equation at the feed side of the membrane reads:

and for the permeate side:

where D is the diffusion coefficient in the membrane (the term D∂c/∂y formulates the diffusive flux along y axis). The CFD model developed contains specific equations for the cases studied, but it can be changed to cover any combination of variations in wall flux, rejection, viscosity and diffusivity.

3.2.1 Summary

CFD methods are valuable tools for process engineers in the design and simulation fields. The differential models, commonly referred to as CFD, aim to solve the Navier–Stokes equation in its discretized form. The advantage of this approach is the ability to use a variety of discretization schemes that suits the geometrical and process needs. The discretization which is dependent on the mesh construction is always a source of error which can be minimized by proper tuning the cell sizes, shapes and types. This formulation of the fluid dynamics approach is not able to predict the thermodynamics of the simulated entity, however. On the other hand, the lattice Boltzmann method, by proper use of collision operator, is able to predict specific thermodynamic behavior of the fluids simulated, e.g. the existence of multiphase system. Such possibility, in order to be incorporated into traditional discrete approach, needs to be accounted for by separate models, for example volume of fluid model of multiphase flow, which further enlarges the formulation and computing demand. One of the key benefits of lattice Boltzmann method for simulating fluid flows as compared to traditional numerical methods is their aptitude to efficiently model interfaces between two or more fluids. Both the fluid dynamics approaches, however totally different, can successfully be incorporated into modeling and simulations of membrane systems.