• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    A review on hollow fiber membrane module towards high separation efficiency: Process modeling in fouling perspective

    2022-09-16 05:23:48XianhuiLiMohammadYounasMashallahRzakazmiQuangVitLyJianxinLi
    Chinese Chemical Letters 2022年8期

    Xianhui Li, Mohammad Younas, Mashallah Rzakazmi, Quang Vit Ly,Jianxin Li

    a Key Laboratory for City Cluster Environmental Safety and Green Development of the Ministry of Education, Institute of Environmental and Ecological Engineering, Guangdong University of Technology, Guangzhou 510006, China

    b Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou), Guangzhou 511458, China

    c State Key Laboratory of Separation Membranes and Membrane Processes, National Center for International Joint Research on Membrane Science and Technology, Tiangong University, Tianjin 300387, China

    d Department of Chemical Engineering, University of Engineering and Technology, P.O.Box 814, Peshawar, Pakistan

    e Faculty of Chemical and Materials Engineering, Shahrood University of Technology, P.O.Box 3619995161, Shahrood, Iran

    f Institute of Research and Development, Duy Tan University, Da nang 550000, Vietnam

    ABSTRACT Hollow fiber microfiltration (MF) and ultrafiltration (UF) membrane processes have been extensively used in water purification and biotechnology.However, complicated filtration hydrodynamics wield a negative influence on fouling mitigation and stability of hollow fiber MF/UF membrane processes.Thus, establishing a mathematical model to understand the membrane processes is essential to guide the optimization of module configurations and to alleviate membrane fouling.Here, we present a comprehensive overview of the hollow fiber MF/UF membrane filtration models developed from different theories.The existing models primarily focus on membrane fouling but rarely on the interactions between the membrane fouling and local filtration hydrodynamics.Therefore, more simplified conceptual models and integrated reduced models need to be built to represent the real filtration behaviors of hollow fiber membranes.Future analyses considering practical requirements including complicated local hydrodynamics and nonuniform membrane properties are suggested to meet the accurate prediction of membrane filtration performance in practical application.This review will inspire the development of high-efficiency hollow fiber membrane modules.

    Keywords:Hollow fiber microfiltration and ultrafiltration process Mathematical modeling Filtration hydrodynamics Membrane fouling Optimization

    1.Introduction

    Membrane technology has been widely applied in various industries such as medical, chemical, energy, environmental and food industries around the world because it only requires less energy during operation and more importantly it is easier to scale up [1,2].Commonly, membrane filtration system is a preferred unit operation in industrial applications owing to its mild operating condition.As the core component of membrane system, the hollow fiber membrane (HFM) module is more preferable over other membrane configurations such as spiral wound and flat sheet,etc.The most attractive advantage is its self-supporting property, leading to a larger effective surface area and higher packing density [3,4].However, compared to the other membrane configurations, the decrease of local flux and trans-membrane pressure (TMP) along flow direction are more prominent in the hollow fiber configuration because of the serious internal pressure drop in the narrow fiber lumen [5,6].Unfortunately, this non-uniform local flux distribution will result in more complicated local hydrodynamics and severer membrane fouling [7,8].It should be noted that the instability of local filtration behavior will be intensified as increasing the membrane module size due to the enhanced non-distribution of local flux distribution.Therefore, much attention has been paid to the optimization of the HFM module in order to realize high efficiency in practical applications [9–11].

    Optimizing operating conditions is a road to enhance the performance of the HFM system, however, this is a time-consuming and expensive job.Mathematical modeling has significant advantages over experimental measurements such as cost-effectiveness,safety, and flexibility [12].In addition, not fully understanding of the transport mechanism occurring inside the HFM system prob-ably overlooks the interaction of important parameters.Therefore,to facilitate a comprehensive consideration during the HFM module design, modeling the membrane separation characteristics is indispensable.However, membrane filtration modes (cross-flow,dead-end, and submerged filtrations) and membrane module configurations (vertical and horizontal configurations) significantly affect the solute transportation [13–15].For instance, in contrast to the dead-end and submerged filtrations in the absence of retentate, the cross-flow filtration enables retentate flowing continuously from the module outlet.This causes a strong shear force on the membrane surface, which benefits fouling mitigation [15].Compared to the horizontal configuration, the vertical configuration suffered from more severe fouling because of more nonuniform local flux distributions [13].The massive difference between these membrane filtration principles results in various theoretical models and numerical approaches.In the last decade, plenty of models have been proposed to predict filtration performance and elucidate fouling mechanisms based on different theories, including mass and momentum transfer, hydrodynamics,etc.[16–18].However, there are few comprehensive reviews regarding the relationship between membrane fouling and local filtration properties of HFM.

    This work aims to synthesize previous studies and current state-of-the-art in modeling studies of hollow fiber microfiltration (MF)/ultrafiltration (UF) membrane processes, which are classified based on their applied theories (Table 1).The filtration phenomenon described by each model will be discussed in the following sections, aiming to guide the selection of the best model for a particular simulation case.Possible improvements of these membrane-based filtration models are discussed and future research scopes are suggested, so as to provide some guidelines concerning the development of appropriate HFM modules for practical applications.

    Table 1 Summary of different models reported in the literature.

    2.Modeling the filtration process

    Fouling is a major and inevitable problem in membrane modules and applications, including HFMs.It reduces permeate flux over time, so as to play its role in reducing cost efficiency [19].Fouling can be defined as the blockage on the pores and surface of the membrane by soluble or suspended matters, such as colloids,biological, inorganic, and organic components [20].Fouling can be solved through chemical, physical, and hydrodynamic cleanings.To fully understand the solute transport mechanisms and interactions between foulants and filtration characteristics, different calculating approaches are proposed to predict the membrane filtration performance.In addition, to accurately extrapolate the filtration of behavior of hollow fiber membranes, the significant parameters of module geometry (e.g.,fiber diameter, fiber length, and number of fibers) have been incorporated into these models.The limitations and improvements of each model are also discussed at the end of each section.

    2.1.Resistance-in-series model

    During filtration of suspension, hydraulic resistance rises significantly owing to the deposition of foulants and the concentration polarization.Darcy’s Law is the basic model that describes the hydraulic permeability,k, for the hollow fiber stating asv=?kP/μ,where μ is the viscosity of the fluid,Pis pressure, andvis the velocity [21].According to the Darcy’s Law, the flux is considered to be controlled by several hydraulic resistances in series (Eq.1)[22,23]:

    whereΔPis the pressure (Pa),Rmis the resistance due to membrane thickness (m?1),Rcis the resistance of the deposited foulants (m?1),Rbis the boundary layer resistance (m?1).

    Generally, four models are used for the description of fouling behaviors depending on the effective particle sizes as well as the membrane pore size as follows:

    (1) The standard blocking model premises an accumulation on the pore walls of the membrane.The membrane permeability is significantly reduced as the pores become constricted.In this case, the decrease of membrane permeability is related to a reduction in the total pore volume.

    (2) The intermediate blocking model premises that a portion of particles block some pores while the rest accumulate on the top of deposited particles.

    (3) The complete blocking model is based on the assumption that the particles have a similar size with the membrane pores,and each particle arriving to the membrane easily blocks pores without superposition of particles.The effective membrane filtration area is thus decreased resulting from the particles deposited on the membrane.

    (4) Cake filtration model premises that particle accumulation results in a permeable cake with increasing thickness.

    These mechanisms have been used separately or integrally to explain experimental results.For instance, Leeet al.proposed a model to describe TMP increase during the filtration of starch solution on the basis of a complete blocking theory in a submerged hollow fiber membrane filtration process [24].It pointed out that fouling easily occurred around the potted area of the membrane close to the outlet of the module, which should be cleaned more frequently than the other parts of the membrane.Another modified model was derived to account for the relationship between the filtration characteristics (feed water characteristics and membrane properties) and fouling mechanisms [25].The simulated results showed that the major fouling mechanism of hydrophobic and hydrophilic membranes for the filtration of raw water was standard blocking and cake formation, respectively.In particular,by separately examining the pressure drop across the cake and membrane, Kinlet al.established a model for prediction of the local filtration behavior along the HFM attached with a compressible cake layer [26].The cake permeability (kc) was expressed as Eq.2.

    whereαcis the specific cake resistance (m?1),pfis the pressure on the feed side (Pa),pois the pressure on the fiber outer wall(Pa),SαandSεare the compressibility index towards the specific cake resistance and porosity, respectively,ρsis the density of suspended particles (kg/m3) andεcis the filtration cake porosity.It demonstrated that the maximum local flux appeared at the place with the highly compressible cake layer.If the resistance is offered by pressure and radial flow in hollow fiber, the total resistance can be calculated by the inverse of hydraulic permeability using the following equation (Eq.3):

    whereLis the fiber length,ΔPis pressure,ro2andri2are the fiber outer and inner radii, respectively, andKis the fiber conductivity[27].

    The concept of resistances-in-series is a convenient way to provide an overview of the fouling contributions in the membrane filtration process.Nevertheless, no details of fouling layers can be given behind the models, since there is a lot of varieties in this calculating approach.For example, current resistances-in-series approaches cannot accurately predict protein filtration behavior because the protein undergoes complicated transformations such as polymerization, aggregation, or coagulation.Therefore, linking fouling to the concentration polarization phenomenon, local hydrodynamics, and force balance of foulants is the significance of accurately predicting real filtration performance and finding the fouling control strategies.

    2.2.Hagen-Poiseuille theory-based model

    The Hagen-Poiseuille equation as a basic hydrodynamic equation is particularly applied to calculate a pressure drop of laminar flow through a rigid pipe [28,29].It can be expressed as Eq.4.

    wherevis the flow velocity (m/s),diis the inner diameter (m),μis the fluid viscosity (Pa.s).

    Therefore, it has been extensively employed to predict the filtration behavior of HFM modules.For example, Carrollet al.simulated the local TMP gradients along the fiber in the double-end and dead-end modules based on Hagen-Poiseuille theory [30,31].They found that the non-uniform distribution of the local flux led to the non-uniform deposition of foulants, but the local flux profile got more uniform with the cake growth.Further, Changet al.developed a mathematical model incorporating the local flux concept into the Hagen-Poiseuille equation to investigate the effect of fiber features (radius, length, and membrane resistance) on filtration resistance and axial pressure drop along the fiber [9].The critical flux is regarded as the flux at which the drag forces of foulants towards the membrane are equal to the shear forces of foulants that transported away from the membrane.It was assumed that as the maximum initial flux at the fiber outlet (Jimax) was lower than the critical flux (Jcr), no deposition occurred and the filtration resistance kept constant throughout the filtration.The local flux distribution was expressed as follows (Eqs.5 and 6) based on Hagen-Poiseuille equation.

    where

    whereLis the fiber length (m),Jmiis the average flux (m/s),Rtotis the overall resistance (m?1).

    Under the condition ofJimax>JcrandJmi

    whereYis the specific productivity of the permeate flow(m3/kWh),φis the fiber packing density,P(I) is the pressure at fiber dead-end (Pa),P(L) is the pressure at fiber outlet (Pa),diis the fiber inner diameter (mm) andδis the fiber thickness (mm).The simulation suggests that the optimal fiber inner radius is 0.20–0.35 mm as the fiber length is 0.5–3.0 m.

    To obtain maximum energy efficiency of submerged hollow fiber membrane module, Yoonet al.developed a model containing the module configuration and packing density as well as the optimum range of the fiber length and diameter for both one- and two-side suctions based on the Hagen-Poiseuille theory [32].In this model, the optimization objective functions were the pressure difference between both ends of the fiber (Df) and the power effi-ciency (Q/P) (Eqs.8 and 9), respectively.

    whereDois the fiber outer diameter (m),εis the area ratio of fiber bundle,dis the distance among fibers (m),Qairis the airflow rate(m3/min),Ais the total footprint (m2) andHis the pressure head(m).

    Results showed that for dead-end filtration mode, the optimum ranges of the fiber length and outer diameter were estimated to be 1.0–1.9 m and 2.2–3.2 mm, respectively, as the operating flux was 20–50 L m?2h?1and the ratio of outer to inner diameters was 2:1[32].For double-end filtration mode, the optimum fiber length and diameter were 1.0–1.9 m and 1.3–2.0 mm under the same operating conditions, respectively.

    Furthermore, a general modeling to predict the flux profile,pressure, and flow rate in a hollow fiber membrane module containing adsorbent beads was developed by Daiet al.[33].The permeate flow rate at any length of the module (Q) was calculated as follows (Eq.10) [33]:

    where the parameters ofA1andA2are the membrane permeation constants, which are related with the local fluxes at the beginning and dead-end of the module, respectively.

    In addition, to reveal the uniformity of TMP profile along the fiber, a parameter called the uniformity factor,β, was incorporated into the model defined as follows (Eq.11):

    where TMP0and TMPLare the TMP at the beginning and end of the membrane module, respectively.Theβrequired is calculated as the maximum value between TMP0and TMPL.A lowerβcorresponds to a more uniform TMP profile.Results showed that this model could successfully predict the TMP and solvent flux profiles[33].

    To investigate the effect of membrane fouling on the local flux,Changet al.further proposed an improved model based on a cake filtration theory [5].The cake filtration theory correlated the cake structure to the electrostatic interaction with a concept of force balance between the inter-particle interactions and the hydraulic forces.The resistance of cake layers (Rcr) formed on the membrane surface was calculated by Eq.12.

    where?is the Happel function,rois the cake solid volume fraction andais the particle radius (m).

    It was found that the foulants deposited on the membrane surface could influence the local flux distribution, which tended to be uniform once the particles were deposited on the membrane[5].In addition, the non-uniform distribution of the cake resistance would intensify as the increase in the operating flux and the decrease in the fiber inner diameter, particle size, and particle zeta potential.Specifically, Wang and Waite paid particular attention to the cake layer growth rate and the local flux distribution on the HFMs using a model developed by combining Darcy law and Hagen-Poiseuille equation [34].It showed that the local flux profile was non-uniform and not constant, and thus resulting in a non-uniform growth rate of foulants, cake layer, and average porosity during HFM filtration.In addition, they also found that the non-uniformity was most significant at the initial filtration but was gradually diminished resulting from cake layer build-up as filtration progressed [34].

    To further investigate the aeration affecting fouling control, Liuet al.[35] simulated the time-dependent filtration performance by introducing an empirical parameter named gas bubbling cleaning efficiency (E) based on the basic method provided by Changet al.[5].Herein,Erepresented the influence of the back-transport force on the foulants away from the membrane surface induced by bubbling, and is related to the nozzle size and the gas flow rate.The cake resistance can be calculated byin the absence of gas bubbling but should be reduced to(1 ?E)under the condition of gas bubbling.The local flux distribution was found to become more uniform with an increasedE[35].AlthoughEpartly reveals the effect of gas bubbling on the local filtration behavior,it is only an empirical parameter.Furthermore, a model for the biological system was created, which incorporated the geometry(e.g., the module layout structure) and hydrodynamics of the system (e.g., the two-phase plug flow in the feed side and hydrodynamics within the fiber) [36].This model showed that high concentration of extracellular polymeric substances strongly increased biofouling, while high concentration of solids enhanced pore blocking and cake layer formation.In addition, although an increase of the aeration rate led to decreasing the thickness of cake layer, it also resulted in an increase of the specific resistance of the cake layer.

    In summary, these models established based on the Hagen-Poiseuille theory provide an important guidance for the module optimization and process control of membrane systems.However,it should be noted that the Hagen-Poiseuille equation is only valid for a horizontal impermeable pipe, which is not the case for membranes that are permeable to water.Hence, the models based on the mass and momentum balances have been developed to solve this issue, elucidating the effect of the transmembrane flow on the TMP and local flux distribution along the permeable HFM[10,16,37].

    2.3.Mass and momentum balances-based model

    To take into account the transmembrane flow to more precisely predict the filtration behavior, a mathematical model was developed by Liet al.based on the momentum balance and complete mass balance [10].A slip parameter (α) was introduced to address the influence of the transmembrane flow on the internal flow resistance (Eqs.13–15),

    where

    Here,fis Fanning frictional factor,RerandReare radial and axial Reynolds number, respectively.αis the dimensionless slip parameter,Kis the membrane permeability (m s?1Pa?1),andvare the axial and radical flow velocity (m/s) inside the fiber, respectively.

    Liet al.reported that by introducingαto address the transmembrane flow, the simulations fit better with the experimental results than the traditional model on the basis of the Hagen-Poiseuille theory [10,37].For instance, the deviation between the experimental and simulated suction pressure calculated by the developed model (1.1%) was less than that (12.8%) obtained by the Hagen-Poiseuille theory as the operating flux was 22.8 L m?2h?1.The non-uniformity of the local flux distribution was found to be positively correlated with the fiber length and operating flux, but negatively correlated with the fiber inner diameter.Particularly, to optimize the fiber diameter and length, the water production effi-ciency (Y) of module as an objective function was calculated (Eq.16) [10]:

    wheregis the gravity constant,Lis the fiber length (m),Lwis the height difference between fiber and air nozzle (m),ρois the feed density (kg/m3),Lhis the distance from the liquid level to the fiber outlet (m),Qis the flow rate (m/s),Qois the aeration rate(m/s),pis the standard atmospheric pressure (Pa),Eis the energy consumption of membrane system (W).

    Fig.1.The influence of fiber diameter and length on the water production efficiency at the different operating flux of (a) 10, (b) 15, (c) 20 and (d) 25 L m?2 h?1, respectively.Copied with permission [10].Copyright 2015, American Institute of Chemical Engineers.

    Fig.2.(a) Schematic of forces excreted on a deposited particle and (b) comparison of the mass and average specific resistance of the cake between the experimental and calculated results.Copied with permission [40].Copyright 2010, WILEY-VCH Verlag GmbH and Co.KGaA, Weinheim.

    As shown in Fig.1a, under the operating flux of 10 L m?2h?1,Yincreased firstly and then reached a plateau with increasing the fiber length and decreasing the fiber inside diameter [10].As the operating flux increased to 15 L m?2h?1(Fig.1b),Yfirstly increased and then reduced with the increase of the fiber length.This trend was more obvious at the operating fluxes of 20 and 25 L m?2h?1(Figs.1c and d).Consequently, selecting a large fiber diameter and a reasonably long fiber length would benefitYat a certain averaged operating flux.

    In addition, to consider the effect of the fluid motion and convention on the permeate flux, Yeh used a complete momentum balance model integrated with the exponential model to simulate the permeate flux [16].The simulated results demonstrated that the operating pressure had a more significant effect than the feed concentration on the permeate flux.The completed mass and momentum balances-based models are lack of foulants properties influencing on filtration performance.The force balance models as other powerful tools have been used to analyze the behaviors of foulant deposition on membrane surface [38–43].

    2.4.Force balance model

    A force-balance model on the basis of the hydrodynamics theories was proposed to study the effects of operating flux, bubble size, and aeration rate on the specific filtration resistance and cake mass in the membrane filtration process [38–40].This model includes the inertial lift force (Fl), the drag forces generated by the permeate flow (Fp), the bubble flow (Fb), the suspension flow (Ft),the net gravitational force (Fg), and the net interparticle force (Fi)(Fig.2a).For example, the increased trend of cake mass (wc) and specific resistance (aav) were predicted by the model, with the theoretical curves very close to the experimental trends (Fig.2b).A narrower particle size distribution in the cake layer with a higher specific resistance was found at an increased aeration rate and a decreased bubble size.These predicted phenomena in accordance with the observations in experiments suggested the accuracy of the model to describe the cake properties.

    To investigate the effect of the particle with a wide size distribution on cake properties and filtration behavior, Lu and Hwang[41] proposed a force balance model for the cross-flow microfiltration (MF).The probabilities of particle deposition and the particle size distribution could be predicted from the theoretical critical frictional angle on the basis of the force analysis.It was assumed that when the frictional angle between particles was larger than the critical angle, the particle would be easily swept away.Otherwise, the particle would deposit stably.The predicted values of the specific resistance and cake porosity agreed very well with the experimental results during filtration of CaCO3solution.Hwanget al.[42] further investigated the effect of particle sizes from submicron to micron on the cross-flow filtration performance.They showed that the size distribution of particle in the cake layer became relatively uniform with increasing the cross-flow velocity, but became non-uniform when increasing the filtration pressure.Meanwhile, a higher filtration pressure resulted in a reduction in the cake porosity but an increase in the cake resistance.An opposite trend of cake properties was observed with the increasing cross-flow velocity.Moreover, to predict the filtration behavior under the aeration condition, a more common case for the submerged HFM system in practical applications, another complete filtration model was developed by combining the force balance model with the gas-liquid two-phase hydrodynamic model [43].It demonstrated that the aeration could aggravate the non-uniform profile of the local flux at the initial stage.However, since a larger local flux would result in more severe fouling, the local flux distribution would be selfadjusted by fouling as the filtration processed to become more uniform.

    Although these models provide deep insights into the filtration characteristics of the HFM, some challenges should be faced to simulate the real filtration process because of the complicated local hydrodynamics especially in the submerged HFM filtration process[44].Since the aeration not only disturbs the local gas-liquid-solid three-phase flow but also induces the movement of the flexible fibers, which intensifies the complication of local hydrodynamics[45], development of a simplified model incorporating the movement of fiber is the significant importance of prediction of real HFM filtration behavior.The local hydrodynamic around fiber not only determines the foulants removal but also affects local flux distribution along fiber [37].Establishment of a reliable hydrodynamics model helps us to comprehensively understand the local hydrodynamic environment around hollow fiber, and would be discussed in following section.

    2.5.Convection and diffusion model

    The convection and diffusion equation is an attractive tool to describe the transport phenomenon occurring during the separation of solutes particularly in the convention-enhanced filtration process [46,47].For instance, Wüpperet al.[48] developed the convection and diffusion models to simulate the local filtration performance of hollow fiber dialyzers.It demonstrated that the transport of low and high molar mass solutes across the high-flux HFM was influenced by diffusion mechanism and diffusion-convection mechanism, respectively.Similarly, Curcioet al.[49] developed a semi-empirical model to differentiate the contributions of diffusion and convection to the transport of different solutes (i.e., the bovine serum albumin, glucose, apo-transferrin, and ammonia).The results showed that for the polysulfone (PSf) HFM, the convection only contributed significantly to the transport of low molecule weight solutes.However, for the polyethersulphone (PES) HFM, the convective contribution increased with the increase of the solute molecule weight.

    To overcome the limitations of classic models, including the concentration polarization model and the resistance-in-series model, a two-dimensional (2D) model on the basis of the convection-diffusion equation was developed to investigate the transport phenomena in UF [50].The predicted permeate fluxes were in excellent agreement with experimental data.The simulated results indicated that the fiber length had a significant influence on the profiles of the local concentration and permeate flux.Meanwhile, the convection-diffusion equations based 2D transport model could help in-deep analysis of the effect of module geometry and membrane transport properties on the efficiency of solute transport [51].The simulations indicated that the enhancement of convective and diffusive solute transport is beneficial for the removal of small-sized molecules by designing module geometry.

    In addition, Mortiet al.[52] provided a quantitative comparison of the convective and diffusive transport properties of homogeneous and asymmetric hemodialysis membranes using transport models.The local solute flux (Ji) of the HFM owning to both convection and diffusion could be expressed as follows (Eq.17):

    whereεbis local porosity,D∞is the Brownian diffusion coeffi-cient (m2/s),Cis the solute concentration (kg/m3), andKDandKCare the hindrance factors for diffusion and convection, respectively,representing the hydrodynamic interaction between the pore wall and solute.

    Fig.3.Variation of pressure along the fiber at the permeate side and at the retentate side for outside-in filtration mode.Copied with permission [18].Copyright 2009, Elsevier.

    The results presented that the clearance of the asymmetric membrane was more dependent on the molecular weight of solute in comparison with a homogeneous membrane [52].The sieving coefficient of the asymmetric membrane was less than those obtained from the homogeneous membrane.Similarly, Bakhshayeshiet al.[53] proposed a theoretical model based on the convection and diffusion models to illustrate the axial variations in TMP and bulk mass transfer coefficient.It confirmed that the TMP profiles and the bulk mass transfer coefficients were related to the flowrate through the fiber lumens and the growth in the concentration polarization boundary layer, respectively.On the basis of this work,a theoretical framework was further developed to investigate the influences of operating conditions (feed and permeate flowrates)and the module geometry (fiber length and number of fibers) on the dextran retention profile in HFM module [54].It demonstrated that the use of a module with a shorter fiber and higher packing density could obtain the best trade-off between backfiltration and concentration polarization.In addition, the profile of dextran retention was insensitive to variations in the flow rates of feed and permeate at the optimal conditions.

    Several studies have used the computational fluid dynamics(CFD) tools to study the influences of operating parameters on the performance for different applications using different geometries of HFM modules [55,56].Usually, finite element approach and finite volume approach are used in CFD simulations.Borsi and Lorain [57] presented a quasi-3D hydrodynamic model coupled with transport, adsorption and attachment equations.The parameters in this model depended directly upon the module characteristics (the number of fibers, membrane thickness,etc.).In addition, this simple model can be used for long-term simulations and allows quick testing and comparison of the filtration behaviors under different filtration parameters.For example, it can easily predict whether the plant could work at a higher operating flux with an improved or at least unchanged hydraulic efficiency.Similarly, the hydrodynamic behavior of the submerged HFM bioreactor at different flow directions and viscosities was explored by Leslieet al.using a porous media model and MBR CFD model [58,59].

    Specifically, Güntheret al.[18] investigated the influence of the packing density on the permeate flow velocity along the fiber using a CFD based model.The results showed that increasing the packing density would result in a non-uniform permeate flux distribution.Especially, when the packing density was very high, filtration flux would decrease dramatically.In contrast, when the packing density was low, the permeate flux was higher at the upper part of the fiber (Fig.3).The influence of the packing density on the spatial deposition of particles along the fiber was further predicted by a numerical model [60], which was in accordance with the trends with respect to the permeate flux.With the highest packing density, particles deposited more easily at the bottom of the fiber for both the inside/out and outside/in filtration modes.This trend was inversed for the case of the lowest packing density.Consequently, a moderate packing density could result in a higher permeate flux and a more homogeneous flux profile along the fiber.

    In addition to the convection-diffusion equation, the concentration polarization model is a powerful and simplified tool to describe the transportation characterization inside the membrane filtration process.

    2.6.Concentration polarization model

    To predict the filtration performance of macromolecules such as protein and dextran in the cross-flow HFM modules, Bakhshayeshiet al.established a mathematical model to simulate the permeate flux on the basis of the concentration polarization model [61].The permeate flux (J) can be calculated by the following expression (Eq.18):

    whereCbis the bulk concentration (mol),Csis the solute concentration on the membrane surface (mol),kmis the mass transfer coefficient,Sois the solute sieving coefficient.

    The fouling layer was assumed as a constant viscosity with a linear fluid velocity profile.This model could successfully predict the permeate fluxes at different particulate suspensions [61].To take into account the effect of the finite cake layer thickness, a modification of the concentration polarization model with shearinduced diffusion was proposed to describe the cake thickness and local permeate flux [62].This model confirmed that a stagnant gel layer in the UF significantly affected the steady permeate flux.A similar model was developed to predict the profile of the cake thickness and filtration flux along the fiber in the cross-flow module by considering the convective flow and the hydrodynamic diffusion of the particle layer [63].It was found that the cake volume fraction and the local TMP varied from the inlet to the outlet of the module.However, the application of these developed models is limited to cases where the permeate flux is independent of TMP.

    To improve the accuracy of prediction, Mondalet al.further combined the concentration polarization model and serial resistance to investigate filtration behaviors in the steady-state laminarflow regime [64].This model was used to determine the transition from the mass- to pressure-transfer-controlled regimes along the fiber.The predicted fluxes were in accordance with the measured values.In addition, by considering the solute diffusivity near the membrane surface, Yehl and Zydney proposed a modified concentration polarization model to predict the permeate flux, which was much closer to the experimental data [65].Coupling the concentration polarization model with local hydrodynamic behavior, such as Navier-Stokes and continuity equation, could provide insight into the relationship between the local fouling and filtration features[66–68].

    2.7.Navier-Stokes and continuity equation

    Flow hydrodynamics through porous media and inside the porous tube can be simulated by the well-known Navier-Stokes and continuity equations, which also inspire researchers to develop models for membrane systems [69,70].For example, a model describing the radial and axial TMP in a cylindrical cross-flow filtration module was proposed on the basis of the continuity and momentum equations [66].In addition to the TMP drop, this model coupled with the concept of shear-induced diffusion and the resistance model could simulate the average permeate fluxwhich can be calculated by using following equation (Eq.19):

    Fig.4.Comparison of the constant flux and pressure modes in terms of the crossflow velocity ratio RU and permeate velocity ratio RV.Copied with permission [71].Copyright 2011, American Institute of Chemical Engineers.

    whereLis the filter length (m),ais the particle radius (m),τwois the shear force (Pa),φois the volume fraction of bulk suspension particle,ηis the relative viscosity (μ(φ)/μo) andis the dimensionless excess particle flux.

    Mondor and Moresoli [66] found that the radial TMP drop decreased from the inlet to the outlet.The difference of TMP drop between the fiber outlet and inlet is related to the fraction of the fluid flowing across the membrane.Further, they provided an extension of the model to the case as a stagnant layer was formed on the hollow fiber membrane surface [67].This model was in consideration of the axial TMP variation.The effect of the stagnant particle layer on the pressure drop across both the fouling layer and the membrane was investigated using this model.It demonstrated that the thickness of stagnant layer increased from the inlet of the fiber, and gradually reached a maximum and then decreased at the outlet of the fiber.In addition, the effect of the fiber radius and particle size on the average permeate flux was further investigated using the developed model [68].The predicted results indicated that the position where a cake began to form gradually moved to the outlet as a result of a reduction of shear-induced diffusion back-transport.In addition, the average permeate flux increased with an increase of lumen radius.These predicted results indicated that both hydrodynamic and back-transport mechanisms simultaneously play a vital role in the cross-flow MF operations[68].

    Furthermore, Kimet al.[71] proposed a perturbation approach of the continuity and Navier-Stokes equations to study the influence of the hydraulic resistance and membrane length on filtration behavior inside the fiber.It was found that the flow velocity inside the fiber exhibited an exponential variation along the fiber.As demonstrated in Fig.4, the constant flux and pressure modes presented a similar velocity profile for the short membrane with the dimensionless transport parameter (βξ)<1.Whereas, the crossflow velocity ratioRUwas always larger than the permeate velocity ratioRVasβξexceeded 1.

    In order to compare the filtration performances of different module configurations, a porous medium model was developed to describe the hydrodynamics of the HFM module operated in a variety of flow configurations, including the dead-end, the closed-shell,and the cross-flow modes as well as the co-current and countercurrent contacting modes [17,72].The fluid continuity and Darcy law were combined to present a set of 2D differential equations.Specifically, the effect of fiber expansion under wet conditions was also combined into the model.Simulations showed that the filtration performance in a cell-free extracapillary space (ECS) was primarily controlled by the membrane resistance.Compared to the cell-packed conditions, the effect of ECS hydraulic permeability on the permeate flux was significant.This mode was a feasible alternative to the counter-current operation.Moreover, Polyakov and Kazenin [73] derived analytical expressions to investigate the effects of flow rate, TMP, and the fiber geometry on the reversible membrane fouling.It was found that both the rectangular and radial depth membrane filtration modules were much more efficient than the existing cross-flow and dead-end devices.Kostoglou and Karabelas [74] focused on the filtration flow in a single fiber and the unidirectional flow in the multi-fiber module.To solve a series of simpler sub-problems, the geometrical and physical simplifications were introduced.The simplification of geometry was usually developed on the basis of the unit cell concept with stochastic parameters.It was found that the decrease of pressure and the increase of flux led to a non-uniform flux distribution.The increase in the ratio of outer/inner diameter would result in an increase of pressure drop in the fiber lumen [74].Another Reynolds averaged Navier-Stokes turbulence modeling approach was incorporated into CFD modeling to simulate the wall shear stresses in two different membrane module configurations [75].It is proven that a modified module configuration with tangential inlet and outlet has better shear stress distribution than the normal module configuration.

    Further, Zhuanget al.[76] systematically studied the flux profile in an outside-in dead-end HFM module using CFD simulation.Three dimensionless numbers were defined to determine the position with the maximum local flux and to quantify the energy utilization and flux distribution.The simulated results indicated that the better uniformity of the local flux profile but higher energy consumption could be obtained by using wide, short, low permeable, and sparsely packed fibers [77].The impact of the membranemodule shell manifold on the initial flow distribution and fouling deposition in an industrial-scale dead-end HFM module was further investigated by Zhuanget al.[78,79] using the CFD model.The flow impingement and wall flow result in the flux mal-distribution,which intensifies the energy consumption.To realize efficient filtration performance, even distributed and large-area inlet holes are suggested to achieve uniform local flow distribution and decrease energy consumption.In addition, the effect of fiber arrangement on the distributions of velocity field was investigated by Wuet al.[80].The staggered arrangement exhibited greater shear stress than the aligned arrangement at the same operating conditions.Meanwhile, the number and arrangement of membrane modules in a turbulence intensified unit can be optimized to achieve a stronger turbulence kinetic energy by using CFD simulations [81].

    Although Navier-Stokes equation based modeling studies provide an accurate prediction of the hydrodynamics and the mass transfer in the membrane filtration process, accurate prediction of local filtration flux may be achieved with the help of the incorporation of different fouling mechanisms.To properly characterize the real hydrodynamic conditions, it is encouraged to integrate multiple models rather than establish a complicated and time-consuming model.

    3.Perspectives and challenges

    Numerous models including empirical and rational calculations have been developed to predict the effect of membrane type, module configurations and operating conditions on the fouling behavior.However, there is still a lack of knowledge on the feasibility of different fouling models and on the fouling mechanisms behind the fouling dynamic.In particular, for some complicated filtration processes (e.g., MBR includes both biological and physical processes), individual model is difficult to simulate their real filtration performance.For instance, classical resistance-in-series theory has shown its potential for fouling modeling, but this work does not consider the mechanisms such as development of fouling layer, activity and deformation of biofilm.Therefore, the integrated models by taking into account the interaction among different parts of the system are highly encouraged.The resistancein-series model could be integrated with biological kinetic models to describe the filtration performance of submerged hollow fiber membrane in MBR [82].The CFD modeling coupled with Hagen-Poiseuille equation provides a quantitative tool to investigate the influences of backwash intensity and duration, and fiber length on effectiveness of backwashing in MBR system [83,84].However, establishment of reliable models to predict the practical large-scale membrane process is still a significant challenge.Further works are therefore recommended to fully understand the behavior of the major components responsible for membrane fouling, which can upgrade the accuracy of filtration modeling predictions.

    Even if satisfactory results have been obtained in short-term operation, the use of long-term running data for validation of those models should be made, which can improve their accuracy in practical application.The modeling task would become more complex when producing a complicated integrated model, which consumes high computational cost.The complexity of integrated model significantly depends on the structure of the selected single models.More efforts should be devoted to the balance of the selected simple models, and the collection of sufficient and appropriate data for model calibration and validation.

    On the one hand, it is the lack of experimental data appropriate for validating the developed fouling model, because the physical phenomena they predict are beyond the scope of the existing monitoring techniques.Meanwhile, some complex hydrodynamic conditions in the real applications are usually simplified in models,e.g., the non-uniformly distributed air bubbles between fibers[85,86], the axial and transversal liquid turbulent flow [87] as well as fiber movement caused by the injected air [88].The validation of these models is therefore needed by means of novelin-situmonitoring techniques (such as nuclear magnetic resonance, ultrasonic time-domain reflectometry) to be developed [89].

    Another study worthwhile considering is the development of cost model based on the basic filtration models discussed in this study.It would help optimize the configuration of membrane system to satisfy the sustainable operation.Like all models, cost model should be operated at dynamic conditions and need to be validated as well.To enhance the applicability, the simulated cost results are suggested to be compared with full-scale running data.Meanwhile, a rigorous calibration procedure is necessary for the application of cost model.

    4.Conclusions

    This work presented an overview and commentary of published modeling studies on MF/UF HFM processes.These studies give us a comprehensive understanding of the relationship between the local filtration performance and membrane fouling in HFM modules, which offer a strategy to improve the membrane module design with enhanced hydrodynamics.However, challenges still exist in practical applications.In general, the empirical model requires a large amount of data for model establishment and has limited predictive capability.Although rational models provide more realistic analysis, they often fail to include all relevant variables under inexact premises.These unavailable data thus make model application difficult or even impractical.Meanwhile, the local filtration behavior in the dynamic membrane filtration process is both complicated and numerous.Therefore, it may not be practical to incorporate all the complexities into a single model.To ensure the accuracy of prediction, future research may be invested in developing a reduced model for data interpretation and process design.Novel designs that enhance filtration uniformity are crucial for sustainable membrane filtration, especially for emerging highlypermeable membranes, as they suffer from more severe membrane fouling and concentration polymerization.In conclusion, for the development of fouling control strategies and appropriate module design, a deep understanding of local filtration behavior is a crucial issue.

    Declaration of competing interest

    The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

    Acknowledgments

    The study is supported by Program for Guangdong Introducing Innovative and Entrepreneurial Teams (No.2019ZT08L213),National Key Research and Development Program of China (No.2020YFA0211003), Key Special Project for Introduced Talents Team of Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou) (No.GML2019ZD0403), and National Natural Science Foundation of China (No.21878230).

    成人午夜高清在线视频| 婷婷色综合大香蕉| av福利片在线观看| 成人欧美大片| 观看免费一级毛片| 欧美丝袜亚洲另类| 精品一区二区三区视频在线| 亚洲国产精品成人综合色| 成人一区二区视频在线观看| 久久精品国产亚洲av涩爱| 欧美精品一区二区大全| 国产亚洲91精品色在线| 在线 av 中文字幕| 老师上课跳d突然被开到最大视频| 日本熟妇午夜| 国产精品福利在线免费观看| 少妇裸体淫交视频免费看高清| 国产精品美女特级片免费视频播放器| 国产成人a区在线观看| 国产白丝娇喘喷水9色精品| 午夜福利网站1000一区二区三区| 国产亚洲91精品色在线| 国产老妇伦熟女老妇高清| 亚洲国产精品国产精品| 久久久久久久久久久免费av| 菩萨蛮人人尽说江南好唐韦庄| 观看美女的网站| 免费黄色在线免费观看| 久久久欧美国产精品| 国产免费一级a男人的天堂| 美女国产视频在线观看| 午夜福利视频1000在线观看| 乱码一卡2卡4卡精品| 久久久久久久久久久免费av| 精品久久久久久久久亚洲| 乱人视频在线观看| 在线免费观看不下载黄p国产| 夜夜爽夜夜爽视频| 国产精品.久久久| 寂寞人妻少妇视频99o| 天天一区二区日本电影三级| 亚洲性久久影院| 国产男女超爽视频在线观看| av网站免费在线观看视频 | 赤兔流量卡办理| 别揉我奶头 嗯啊视频| 亚洲精品一二三| 国产午夜福利久久久久久| 久久久久久九九精品二区国产| 亚洲av.av天堂| 亚洲精品乱久久久久久| 精品国产三级普通话版| 精品国产露脸久久av麻豆 | 日本色播在线视频| 日本色播在线视频| 男女下面进入的视频免费午夜| 久久99精品国语久久久| 我的老师免费观看完整版| 午夜福利高清视频| 午夜福利在线观看免费完整高清在| 九草在线视频观看| 国产伦精品一区二区三区四那| 亚洲图色成人| 99久久人妻综合| 十八禁国产超污无遮挡网站| 亚洲性久久影院| 白带黄色成豆腐渣| 成人毛片60女人毛片免费| 午夜福利高清视频| www.色视频.com| 精品午夜福利在线看| 777米奇影视久久| 一级毛片电影观看| 国产精品久久久久久av不卡| 午夜激情欧美在线| 国产乱人偷精品视频| 国产 一区 欧美 日韩| 成人综合一区亚洲| 三级经典国产精品| 亚洲精品aⅴ在线观看| 22中文网久久字幕| 亚洲欧美日韩无卡精品| 又粗又硬又长又爽又黄的视频| 国国产精品蜜臀av免费| 欧美日韩视频高清一区二区三区二| 久久精品久久精品一区二区三区| 午夜福利视频1000在线观看| 欧美日韩亚洲高清精品| 亚洲国产日韩欧美精品在线观看| 日本午夜av视频| 国产黄片视频在线免费观看| 欧美成人午夜免费资源| 成年版毛片免费区| 男人舔女人下体高潮全视频| 亚洲va在线va天堂va国产| 日日啪夜夜爽| 永久网站在线| 日韩强制内射视频| 欧美成人精品欧美一级黄| 国产成人freesex在线| 国产精品无大码| 听说在线观看完整版免费高清| 亚洲人与动物交配视频| 天堂网av新在线| 精品国产一区二区三区久久久樱花 | 女人十人毛片免费观看3o分钟| 最近2019中文字幕mv第一页| 久99久视频精品免费| 乱系列少妇在线播放| 欧美高清性xxxxhd video| 精品一区二区三卡| 国产成人一区二区在线| 亚洲精品一二三| 国产成人精品福利久久| 黄色日韩在线| 国产有黄有色有爽视频| 国产免费福利视频在线观看| 熟女电影av网| 国产成人精品久久久久久| 亚洲真实伦在线观看| 男人舔女人下体高潮全视频| 久久久久久久久大av| 亚洲av.av天堂| 乱人视频在线观看| 春色校园在线视频观看| 免费看不卡的av| 欧美日韩在线观看h| 久久精品久久精品一区二区三区| 色尼玛亚洲综合影院| 大陆偷拍与自拍| 美女黄网站色视频| 秋霞在线观看毛片| 成年版毛片免费区| 国国产精品蜜臀av免费| 免费黄网站久久成人精品| 国产精品一区二区在线观看99 | 国产精品国产三级国产专区5o| 伊人久久精品亚洲午夜| 在线观看一区二区三区| 国模一区二区三区四区视频| 国产一区有黄有色的免费视频 | 国产精品久久视频播放| 国产欧美日韩精品一区二区| 国产黄频视频在线观看| 国产黄片视频在线免费观看| 欧美成人午夜免费资源| 男女啪啪激烈高潮av片| 日韩av在线大香蕉| 欧美97在线视频| 精品人妻偷拍中文字幕| 日日摸夜夜添夜夜爱| 色尼玛亚洲综合影院| 日本黄大片高清| 国国产精品蜜臀av免费| 国产久久久一区二区三区| 九九在线视频观看精品| 九色成人免费人妻av| 久久人人爽人人片av| 黄色配什么色好看| 天堂√8在线中文| 在线免费观看不下载黄p国产| 十八禁网站网址无遮挡 | av在线老鸭窝| 国产精品综合久久久久久久免费| 黑人高潮一二区| 超碰av人人做人人爽久久| 亚洲av免费高清在线观看| 尤物成人国产欧美一区二区三区| 伊人久久精品亚洲午夜| 国产av码专区亚洲av| av国产免费在线观看| 激情 狠狠 欧美| 国产视频首页在线观看| 亚洲最大成人av| 床上黄色一级片| 嫩草影院精品99| 性插视频无遮挡在线免费观看| 婷婷色综合www| 国产人妻一区二区三区在| 国产视频内射| 男女那种视频在线观看| 亚洲人成网站在线观看播放| 国产精品一区www在线观看| 欧美成人a在线观看| 韩国高清视频一区二区三区| 国产色婷婷99| 精品人妻视频免费看| 国产精品久久久久久av不卡| 亚洲精品一二三| 亚洲美女视频黄频| 久久精品夜色国产| 国产精品1区2区在线观看.| 久久人人爽人人爽人人片va| 18禁动态无遮挡网站| 欧美三级亚洲精品| 观看美女的网站| 免费黄频网站在线观看国产| 男女边摸边吃奶| 卡戴珊不雅视频在线播放| 亚洲av成人精品一区久久| 国产精品福利在线免费观看| 男人爽女人下面视频在线观看| 搞女人的毛片| 久久久国产一区二区| 久久草成人影院| 国产一级毛片七仙女欲春2| 欧美一区二区亚洲| 国产综合懂色| 国产真实伦视频高清在线观看| 18禁在线无遮挡免费观看视频| 日韩视频在线欧美| 最新中文字幕久久久久| 在线免费观看不下载黄p国产| 免费无遮挡裸体视频| xxx大片免费视频| 亚洲精品中文字幕在线视频 | 国产视频首页在线观看| av在线老鸭窝| 热99在线观看视频| 女人被狂操c到高潮| 全区人妻精品视频| 女人十人毛片免费观看3o分钟| 国产成人a区在线观看| 国产一级毛片七仙女欲春2| 欧美日韩视频高清一区二区三区二| 亚洲天堂国产精品一区在线| 91aial.com中文字幕在线观看| 成人无遮挡网站| 国产在线男女| 国产精品国产三级国产av玫瑰| 97精品久久久久久久久久精品| 天天一区二区日本电影三级| 久久久精品94久久精品| 亚洲久久久久久中文字幕| ponron亚洲| 亚洲欧美日韩卡通动漫| 国产伦理片在线播放av一区| 国产成人精品婷婷| 亚洲av不卡在线观看| 精品一区二区三区人妻视频| 在线观看免费高清a一片| 国产淫语在线视频| 最近最新中文字幕免费大全7| 99热6这里只有精品| 成人亚洲精品一区在线观看 | 中文在线观看免费www的网站| 老女人水多毛片| 亚洲色图av天堂| 国产伦精品一区二区三区视频9| 亚洲婷婷狠狠爱综合网| 99久国产av精品| 国产毛片a区久久久久| 白带黄色成豆腐渣| 日韩人妻高清精品专区| 久久久久精品久久久久真实原创| 校园人妻丝袜中文字幕| 狠狠精品人妻久久久久久综合| 熟妇人妻久久中文字幕3abv| 亚洲精品,欧美精品| 五月伊人婷婷丁香| 一区二区三区四区激情视频| 免费在线观看成人毛片| 99久久精品一区二区三区| 天天一区二区日本电影三级| 日本午夜av视频| 又粗又硬又长又爽又黄的视频| 日日撸夜夜添| 精品一区在线观看国产| 亚洲四区av| 亚洲精品乱码久久久v下载方式| 国产精品麻豆人妻色哟哟久久 | 久久久久久国产a免费观看| 毛片一级片免费看久久久久| 一级黄片播放器| 99热这里只有精品一区| 久久精品国产亚洲av涩爱| av在线观看视频网站免费| 欧美一区二区亚洲| 人人妻人人看人人澡| 国产 亚洲一区二区三区 | 秋霞在线观看毛片| 午夜福利高清视频| 在线免费十八禁| 白带黄色成豆腐渣| 亚洲欧美日韩卡通动漫| 亚洲最大成人中文| 美女xxoo啪啪120秒动态图| 成年女人在线观看亚洲视频 | 大香蕉97超碰在线| 成人亚洲欧美一区二区av| 免费大片18禁| 日韩一区二区视频免费看| www.色视频.com| 欧美日韩亚洲高清精品| 国产 一区精品| 国产精品美女特级片免费视频播放器| 亚洲精品aⅴ在线观看| 最后的刺客免费高清国语| 黄色配什么色好看| 两个人视频免费观看高清| 国产亚洲精品av在线| 亚洲熟妇中文字幕五十中出| 亚洲人成网站高清观看| 亚洲av男天堂| 亚洲图色成人| 99久久精品热视频| 亚洲精品一区蜜桃| 国内揄拍国产精品人妻在线| 国产91av在线免费观看| 免费观看a级毛片全部| 国产人妻一区二区三区在| 欧美日韩一区二区视频在线观看视频在线 | 国产 亚洲一区二区三区 | 欧美日韩亚洲高清精品| 亚洲成人中文字幕在线播放| 99久国产av精品| 国产亚洲5aaaaa淫片| 国产淫片久久久久久久久| 欧美日韩综合久久久久久| 综合色丁香网| 免费看光身美女| 亚洲av电影在线观看一区二区三区 | 亚洲av中文字字幕乱码综合| 欧美精品国产亚洲| 成人综合一区亚洲| 久久久国产一区二区| 国产色爽女视频免费观看| 久久99热这里只有精品18| 国产精品久久久久久久久免| 深夜a级毛片| 国产亚洲av嫩草精品影院| 欧美日本视频| 街头女战士在线观看网站| 三级经典国产精品| 精品久久久久久电影网| 亚洲精品成人av观看孕妇| 精品酒店卫生间| 久久99热这里只有精品18| 一个人免费在线观看电影| 免费观看性生交大片5| 亚洲在线观看片| 热99在线观看视频| 一二三四中文在线观看免费高清| 高清在线视频一区二区三区| 中文天堂在线官网| 少妇裸体淫交视频免费看高清| 韩国高清视频一区二区三区| 国产午夜精品一二区理论片| 高清午夜精品一区二区三区| 亚洲欧洲日产国产| 蜜桃亚洲精品一区二区三区| 高清午夜精品一区二区三区| 美女黄网站色视频| 毛片一级片免费看久久久久| 91精品一卡2卡3卡4卡| 男女啪啪激烈高潮av片| or卡值多少钱| 国产 一区 欧美 日韩| 美女主播在线视频| 日韩 亚洲 欧美在线| 九九在线视频观看精品| 国精品久久久久久国模美| 中文字幕av在线有码专区| 婷婷六月久久综合丁香| 午夜激情欧美在线| 国产白丝娇喘喷水9色精品| 麻豆国产97在线/欧美| 亚洲av电影不卡..在线观看| 永久免费av网站大全| 好男人在线观看高清免费视频| 秋霞伦理黄片| 久久久亚洲精品成人影院| 最近视频中文字幕2019在线8| 久久精品夜色国产| 成人高潮视频无遮挡免费网站| 看非洲黑人一级黄片| 69av精品久久久久久| 国产成人精品一,二区| freevideosex欧美| 视频中文字幕在线观看| 内地一区二区视频在线| 免费无遮挡裸体视频| 国产成人精品婷婷| 亚洲精品中文字幕在线视频 | 80岁老熟妇乱子伦牲交| 国产精品av视频在线免费观看| 国产片特级美女逼逼视频| 国产亚洲精品av在线| 成人国产麻豆网| 插阴视频在线观看视频| 国产日韩欧美在线精品| 寂寞人妻少妇视频99o| 一区二区三区四区激情视频| 午夜精品一区二区三区免费看| 免费高清在线观看视频在线观看| 国产免费视频播放在线视频 | 99热网站在线观看| 秋霞在线观看毛片| 人妻制服诱惑在线中文字幕| 在线免费观看的www视频| 亚洲国产色片| 久久久久国产网址| 97超视频在线观看视频| 亚洲av免费高清在线观看| 在线播放无遮挡| 久久99热这里只有精品18| 亚洲av免费高清在线观看| 成年版毛片免费区| 校园人妻丝袜中文字幕| 成人亚洲精品一区在线观看 | 久久精品夜夜夜夜夜久久蜜豆| 欧美极品一区二区三区四区| 日日撸夜夜添| 免费人成在线观看视频色| 久久久久久久午夜电影| 久久久久久久久久成人| 国产精品久久久久久av不卡| 免费观看在线日韩| 高清欧美精品videossex| 我的老师免费观看完整版| 欧美xxxx黑人xx丫x性爽| 欧美激情国产日韩精品一区| 国产av在哪里看| 精品国产露脸久久av麻豆 | 国产一区二区亚洲精品在线观看| 18禁动态无遮挡网站| 亚洲精品aⅴ在线观看| 久久久色成人| 国产精品久久久久久久电影| 国语对白做爰xxxⅹ性视频网站| 久久久久久久久久人人人人人人| 一二三四中文在线观看免费高清| 日韩,欧美,国产一区二区三区| 成人漫画全彩无遮挡| 久久久欧美国产精品| 狂野欧美激情性xxxx在线观看| 国产欧美日韩精品一区二区| 一级a做视频免费观看| av专区在线播放| 国产成人精品婷婷| 免费观看精品视频网站| 午夜激情久久久久久久| 亚洲三级黄色毛片| 777米奇影视久久| 国产色婷婷99| 日韩av免费高清视频| 成人一区二区视频在线观看| 日本猛色少妇xxxxx猛交久久| 成人美女网站在线观看视频| 免费av毛片视频| 久久久久性生活片| 精品不卡国产一区二区三区| 国产欧美另类精品又又久久亚洲欧美| 国产精品美女特级片免费视频播放器| 欧美最新免费一区二区三区| 免费人成在线观看视频色| 国产高清有码在线观看视频| 成人亚洲精品av一区二区| 波野结衣二区三区在线| 亚洲av免费在线观看| 成人性生交大片免费视频hd| 小蜜桃在线观看免费完整版高清| 3wmmmm亚洲av在线观看| 亚洲va在线va天堂va国产| 欧美日本视频| 丰满乱子伦码专区| 精品一区二区三卡| 亚洲一级一片aⅴ在线观看| 熟妇人妻不卡中文字幕| 97人妻精品一区二区三区麻豆| 美女高潮的动态| 亚洲精品日韩在线中文字幕| 国产人妻一区二区三区在| 国内精品美女久久久久久| 高清欧美精品videossex| 国产老妇女一区| videossex国产| .国产精品久久| 大片免费播放器 马上看| 看免费成人av毛片| 国产综合精华液| 欧美+日韩+精品| 国产精品女同一区二区软件| 欧美日本视频| 日本爱情动作片www.在线观看| 亚洲成人av在线免费| 成年人午夜在线观看视频 | 日本免费a在线| 国产精品一区二区在线观看99 | 国产av码专区亚洲av| 国产麻豆成人av免费视频| 五月伊人婷婷丁香| 夜夜爽夜夜爽视频| 日本av手机在线免费观看| 三级男女做爰猛烈吃奶摸视频| 人妻少妇偷人精品九色| 看非洲黑人一级黄片| 国产成人a区在线观看| av在线天堂中文字幕| 久久久久久久久久久免费av| 午夜爱爱视频在线播放| 男人舔女人下体高潮全视频| 久久精品熟女亚洲av麻豆精品 | 97超碰精品成人国产| 久久6这里有精品| 亚洲电影在线观看av| 国产视频内射| 日韩大片免费观看网站| 久久久久免费精品人妻一区二区| 一区二区三区乱码不卡18| 亚洲av国产av综合av卡| 日日啪夜夜爽| 狂野欧美激情性xxxx在线观看| h日本视频在线播放| 亚洲精品日本国产第一区| 国产 一区精品| 亚洲人成网站高清观看| 天美传媒精品一区二区| 国产精品久久久久久久电影| 欧美丝袜亚洲另类| 国产 亚洲一区二区三区 | 亚洲一级一片aⅴ在线观看| 免费少妇av软件| 成人特级av手机在线观看| 国产午夜精品论理片| 一夜夜www| 国产黄片美女视频| 99久久精品国产国产毛片| 亚洲第一区二区三区不卡| 老司机影院成人| 日韩精品青青久久久久久| 韩国av在线不卡| 中文天堂在线官网| 久久99热6这里只有精品| 热99在线观看视频| 99热这里只有是精品50| 亚洲欧洲日产国产| 国产精品99久久久久久久久| 两个人视频免费观看高清| 99视频精品全部免费 在线| 国产亚洲av嫩草精品影院| 91久久精品电影网| 99热全是精品| 男女下面进入的视频免费午夜| 乱码一卡2卡4卡精品| 777米奇影视久久| 小蜜桃在线观看免费完整版高清| 激情 狠狠 欧美| 麻豆成人av视频| 国产av国产精品国产| 在线a可以看的网站| 99热全是精品| 日韩欧美三级三区| 免费看光身美女| 在线免费十八禁| 热99在线观看视频| 亚洲精华国产精华液的使用体验| 亚洲va在线va天堂va国产| 中文欧美无线码| 免费少妇av软件| 日日摸夜夜添夜夜爱| 日韩人妻高清精品专区| 精品久久久久久久久久久久久| or卡值多少钱| 免费观看在线日韩| 久久久久久久大尺度免费视频| 国产黄a三级三级三级人| 网址你懂的国产日韩在线| 99热全是精品| 亚洲精品色激情综合| 久久精品国产亚洲av涩爱| 久久久久精品性色| 又大又黄又爽视频免费| 精品久久久久久久人妻蜜臀av| 亚洲国产最新在线播放| 亚洲自拍偷在线| 青春草视频在线免费观看| 国产成人a∨麻豆精品| 汤姆久久久久久久影院中文字幕 | 日产精品乱码卡一卡2卡三| 久久久a久久爽久久v久久| 69人妻影院| 久久久亚洲精品成人影院| 国产在视频线在精品| 综合色丁香网| 精品一区二区三区人妻视频| 欧美高清性xxxxhd video| 嫩草影院精品99| 国产一级毛片七仙女欲春2| 中文字幕免费在线视频6| 搡老乐熟女国产| 熟女电影av网| 久久久久久久大尺度免费视频| 成人特级av手机在线观看| 欧美日韩一区二区视频在线观看视频在线 | 最后的刺客免费高清国语| 少妇人妻一区二区三区视频| 国产探花在线观看一区二区| 精品久久久久久久末码| 只有这里有精品99| av黄色大香蕉| 中国国产av一级| 永久免费av网站大全| 高清在线视频一区二区三区| 国产黄片视频在线免费观看| ponron亚洲| 观看免费一级毛片| 69av精品久久久久久| 日产精品乱码卡一卡2卡三| 中文字幕久久专区| 一区二区三区免费毛片| 亚洲乱码一区二区免费版| 久久这里有精品视频免费| 国产在视频线精品| 男女下面进入的视频免费午夜|