• <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| 久久久久九九精品影院| 亚洲成人久久爱视频| 国产在视频线在精品| 一个人看视频在线观看www免费 | 首页视频小说图片口味搜索| 搡老熟女国产l中国老女人| 两个人的视频大全免费| 亚洲精品国产精品久久久不卡| 一个人观看的视频www高清免费观看| 最近在线观看免费完整版| 99久久成人亚洲精品观看| 国产毛片a区久久久久| 亚洲精品一区av在线观看| 亚洲,欧美精品.| 欧美黑人巨大hd| 欧美3d第一页| 99久久99久久久精品蜜桃| 国产高清视频在线播放一区| 成年女人永久免费观看视频| 一进一出抽搐动态| 国产av麻豆久久久久久久| 在线天堂最新版资源| 亚洲成人中文字幕在线播放| 欧美在线一区亚洲| 丰满人妻熟妇乱又伦精品不卡| 国内精品一区二区在线观看| 美女免费视频网站| 丝袜美腿在线中文| 丰满人妻熟妇乱又伦精品不卡| 国产乱人伦免费视频| 91av网一区二区| 听说在线观看完整版免费高清| 麻豆成人午夜福利视频| 麻豆一二三区av精品| 精品福利观看| 啦啦啦免费观看视频1| 国产精品女同一区二区软件 | 亚洲精品久久国产高清桃花| 精品99又大又爽又粗少妇毛片 | 看免费av毛片| 1000部很黄的大片| 国产黄a三级三级三级人| 色精品久久人妻99蜜桃| 亚洲天堂国产精品一区在线| av福利片在线观看| 国产91精品成人一区二区三区| 91久久精品国产一区二区成人 | 精品一区二区三区人妻视频| 亚洲精品粉嫩美女一区| 老司机午夜十八禁免费视频| 午夜免费男女啪啪视频观看 | 亚洲av电影不卡..在线观看| 三级男女做爰猛烈吃奶摸视频| 精品国产三级普通话版| 俺也久久电影网| 黄色丝袜av网址大全| 手机成人av网站| 91av网一区二区| 男女之事视频高清在线观看| 99热这里只有精品一区| 成年女人看的毛片在线观看| 欧美黑人欧美精品刺激| av女优亚洲男人天堂| 亚洲av成人精品一区久久| 国产亚洲精品一区二区www| 欧美区成人在线视频| 熟女电影av网| 在线观看66精品国产| 亚洲专区国产一区二区| 国产精品99久久99久久久不卡| 精品不卡国产一区二区三区| 成人特级黄色片久久久久久久| 精品人妻偷拍中文字幕| 日本与韩国留学比较| 亚洲国产精品999在线| 久久99热这里只有精品18| 亚洲av日韩精品久久久久久密| 国产精品女同一区二区软件 | 国产av在哪里看| 国产精品,欧美在线| 亚洲狠狠婷婷综合久久图片| 老司机午夜福利在线观看视频| 757午夜福利合集在线观看| 欧美在线黄色| 最近最新免费中文字幕在线| 亚洲av成人精品一区久久| 国产免费一级a男人的天堂| 久久精品影院6| 亚洲人成电影免费在线| 午夜日韩欧美国产| 久久国产乱子伦精品免费另类| 超碰av人人做人人爽久久 | 此物有八面人人有两片| 亚洲av第一区精品v没综合| 老司机福利观看| 手机成人av网站| 国产欧美日韩一区二区三| 18禁裸乳无遮挡免费网站照片| 观看免费一级毛片| 国产高潮美女av| 琪琪午夜伦伦电影理论片6080| 免费在线观看亚洲国产| 免费无遮挡裸体视频| 免费av毛片视频| 美女cb高潮喷水在线观看| 午夜精品一区二区三区免费看| 日韩人妻高清精品专区| 免费看美女性在线毛片视频| 日本三级黄在线观看| 欧美区成人在线视频| 免费高清视频大片| 热99re8久久精品国产| 欧美+日韩+精品| 亚洲精品成人久久久久久| 免费人成视频x8x8入口观看| 国产高清三级在线| 欧美又色又爽又黄视频| 天天躁日日操中文字幕| 高潮久久久久久久久久久不卡| 国产精品综合久久久久久久免费| 深爱激情五月婷婷| 一本精品99久久精品77| 两个人的视频大全免费| 免费观看人在逋| 内地一区二区视频在线| 精品一区二区三区人妻视频| 性色av乱码一区二区三区2| 精品99又大又爽又粗少妇毛片 | 亚洲色图av天堂| 国产私拍福利视频在线观看| 黄色片一级片一级黄色片| 91九色精品人成在线观看| 精华霜和精华液先用哪个| 一个人观看的视频www高清免费观看| 最近最新中文字幕大全电影3| 老汉色∧v一级毛片| 久久久国产精品麻豆| 成人鲁丝片一二三区免费| 手机成人av网站| 亚洲av免费在线观看| 99在线视频只有这里精品首页| 真人做人爱边吃奶动态| 国产精品香港三级国产av潘金莲| 99在线视频只有这里精品首页| 久久精品国产综合久久久| 少妇丰满av| 高清毛片免费观看视频网站| 亚洲国产日韩欧美精品在线观看 | 亚洲欧美日韩高清在线视频| 日韩成人在线观看一区二区三区| 免费观看的影片在线观看| 99久久精品一区二区三区| 日本与韩国留学比较| 久久精品91无色码中文字幕| 中文字幕高清在线视频| 久久这里只有精品中国| 国产麻豆成人av免费视频| 欧美精品啪啪一区二区三区| 19禁男女啪啪无遮挡网站| 少妇的丰满在线观看| 久久性视频一级片| 久久久久免费精品人妻一区二区| 国产精品,欧美在线| 99热6这里只有精品| 成人高潮视频无遮挡免费网站| 高清日韩中文字幕在线| 午夜福利视频1000在线观看| 一本综合久久免费| АⅤ资源中文在线天堂| 欧美黄色淫秽网站| 我的老师免费观看完整版| 亚洲狠狠婷婷综合久久图片| 真实男女啪啪啪动态图| 亚洲精品日韩av片在线观看 | 又黄又爽又免费观看的视频| 中文字幕人妻丝袜一区二区| 老司机福利观看| 日本黄大片高清| 亚洲av电影在线进入| 日本三级黄在线观看| 老熟妇仑乱视频hdxx| 一个人免费在线观看电影| 国产 一区 欧美 日韩| 老熟妇仑乱视频hdxx| а√天堂www在线а√下载| 制服人妻中文乱码| 日韩精品青青久久久久久| 欧美激情在线99| 成人18禁在线播放| 精品久久久久久,| 制服丝袜大香蕉在线| 一进一出抽搐动态| 亚洲电影在线观看av| 久久国产精品影院| 国内精品美女久久久久久| 又粗又爽又猛毛片免费看| 国产欧美日韩一区二区三| av在线天堂中文字幕| 亚洲欧美日韩高清专用| 成人18禁在线播放| 国产毛片a区久久久久| 久久久久亚洲av毛片大全| 免费高清视频大片| 午夜免费激情av| 99久国产av精品| 成人精品一区二区免费| 女同久久另类99精品国产91| 美女免费视频网站| 99国产综合亚洲精品| 免费在线观看亚洲国产| 精品久久久久久久人妻蜜臀av| 一区二区三区高清视频在线| 国产精品一区二区三区四区久久| 狂野欧美白嫩少妇大欣赏| 国产毛片a区久久久久| 亚洲精品粉嫩美女一区| 熟女少妇亚洲综合色aaa.| 午夜福利高清视频| 嫩草影院入口| 国产探花极品一区二区| 亚洲内射少妇av| 亚洲男人的天堂狠狠| 久久精品人妻少妇| 美女cb高潮喷水在线观看| 又黄又粗又硬又大视频| 亚洲无线在线观看| 亚洲av成人精品一区久久| 亚洲成人中文字幕在线播放| 中文字幕精品亚洲无线码一区| 国产乱人伦免费视频| 亚洲狠狠婷婷综合久久图片| 国产 一区 欧美 日韩| 在线视频色国产色| 99热6这里只有精品| 久久香蕉精品热| 久久久久久大精品| 脱女人内裤的视频| 男女下面进入的视频免费午夜| 精品一区二区三区视频在线 | 亚洲av免费高清在线观看| 国产精品乱码一区二三区的特点| 淫秽高清视频在线观看| 怎么达到女性高潮| 国产精品免费一区二区三区在线| 国产精品久久久久久久久免 | 一级a爱片免费观看的视频| 精品一区二区三区视频在线 | 操出白浆在线播放| 久久精品91蜜桃| 天美传媒精品一区二区| 欧美日本视频| 色综合站精品国产| 国产探花在线观看一区二区| 国产精品,欧美在线| 中出人妻视频一区二区| 亚洲精品在线观看二区| 1024手机看黄色片| 亚洲电影在线观看av| 99国产精品一区二区蜜桃av| 午夜日韩欧美国产| 琪琪午夜伦伦电影理论片6080| 精品久久久久久久久久久久久| 好男人电影高清在线观看| 亚洲国产欧美人成| 亚洲不卡免费看| 亚洲精品在线观看二区| 国产成人aa在线观看| 网址你懂的国产日韩在线| 丰满的人妻完整版| 亚洲成a人片在线一区二区| 日本成人三级电影网站| 国产欧美日韩精品亚洲av| 亚洲人成伊人成综合网2020| 91字幕亚洲| 青草久久国产| 成人国产综合亚洲| 亚洲欧美一区二区三区黑人| 亚洲国产日韩欧美精品在线观看 | 国产亚洲精品一区二区www| 长腿黑丝高跟| 日本一本二区三区精品| 99久久精品热视频| 国产精品 欧美亚洲| 色在线成人网| 国内精品久久久久久久电影| 日韩免费av在线播放| 午夜老司机福利剧场| 天天一区二区日本电影三级| 国产伦人伦偷精品视频| 久久精品国产综合久久久| 国产av不卡久久| 中文字幕人妻熟人妻熟丝袜美 | 两个人看的免费小视频| 一进一出好大好爽视频| 高清在线国产一区| 久久亚洲精品不卡| 亚洲avbb在线观看| 久久精品国产自在天天线| bbb黄色大片| 亚洲中文字幕一区二区三区有码在线看| 欧美黑人巨大hd| 日韩国内少妇激情av| 国产一区二区三区在线臀色熟女| 久久久久久久久中文| 十八禁网站免费在线| 国产精品99久久久久久久久| 18禁在线播放成人免费| 丁香六月欧美| 亚洲av成人av| 免费观看精品视频网站| 男女床上黄色一级片免费看| 亚洲av成人精品一区久久| 色综合亚洲欧美另类图片| 国产成人欧美在线观看| 此物有八面人人有两片| 在线观看日韩欧美| 91久久精品电影网| 国产亚洲av嫩草精品影院| 亚洲av五月六月丁香网| 夜夜爽天天搞| 国产精品一区二区三区四区免费观看 | 欧美xxxx黑人xx丫x性爽| 欧美+日韩+精品| 中文字幕熟女人妻在线| a级一级毛片免费在线观看| 日韩欧美三级三区| 亚洲在线自拍视频| 岛国视频午夜一区免费看| 丰满人妻熟妇乱又伦精品不卡| 色噜噜av男人的天堂激情| 久久精品国产99精品国产亚洲性色| 伊人久久大香线蕉亚洲五| av女优亚洲男人天堂| 999久久久精品免费观看国产| 成人无遮挡网站| h日本视频在线播放| 欧美zozozo另类| 99国产精品一区二区三区| 一个人免费在线观看的高清视频| 波多野结衣巨乳人妻| 有码 亚洲区| 免费观看精品视频网站| 有码 亚洲区| 啦啦啦观看免费观看视频高清| 宅男免费午夜| 国产aⅴ精品一区二区三区波| 精品国产三级普通话版| www.熟女人妻精品国产| www.色视频.com| 亚洲人成网站在线播| 18禁黄网站禁片免费观看直播| 两人在一起打扑克的视频| 一进一出抽搐gif免费好疼| 丝袜美腿在线中文| 亚洲欧美日韩高清在线视频| 91久久精品电影网| 淫秽高清视频在线观看| 在线播放国产精品三级| 国产精品综合久久久久久久免费| 男女之事视频高清在线观看| 亚洲av电影不卡..在线观看| 国产爱豆传媒在线观看| 亚洲黑人精品在线| 日日摸夜夜添夜夜添小说| 美女大奶头视频| tocl精华| 成人高潮视频无遮挡免费网站| 天堂√8在线中文| 久久精品国产99精品国产亚洲性色| 亚洲精品在线美女| 露出奶头的视频| 久久精品国产99精品国产亚洲性色| 在线观看66精品国产| 国产一区二区在线av高清观看| 成年免费大片在线观看| 久久国产精品影院| 19禁男女啪啪无遮挡网站| 国产成人欧美在线观看| 大型黄色视频在线免费观看| 麻豆成人午夜福利视频| 两个人的视频大全免费| 亚洲欧美日韩高清专用| 老汉色∧v一级毛片| a在线观看视频网站| 亚洲精品一卡2卡三卡4卡5卡| 亚洲黑人精品在线| 国产乱人伦免费视频| 制服丝袜大香蕉在线| 中文资源天堂在线| 99国产综合亚洲精品| 看黄色毛片网站| 亚洲国产欧洲综合997久久,| 天堂影院成人在线观看| 国产蜜桃级精品一区二区三区| 99热这里只有精品一区| 国产成人系列免费观看| 91麻豆av在线| 国产野战对白在线观看| 国产69精品久久久久777片| 午夜福利在线在线| 欧美中文日本在线观看视频| 午夜福利18| 国产成人av教育| 亚洲精品亚洲一区二区| 亚洲人与动物交配视频| 日韩免费av在线播放| 成人性生交大片免费视频hd| 国产精品综合久久久久久久免费| 在线观看免费午夜福利视频| 人人妻,人人澡人人爽秒播| 在线观看日韩欧美| 91字幕亚洲| 一进一出抽搐动态| 97碰自拍视频| 夜夜夜夜夜久久久久| 波多野结衣高清无吗| 国产精品美女特级片免费视频播放器| 欧美一区二区精品小视频在线| 久久九九热精品免费| 精品国产亚洲在线| 欧美黑人巨大hd| 少妇的逼水好多| 有码 亚洲区| a级一级毛片免费在线观看| 久久香蕉国产精品| 一级毛片高清免费大全| 久久久国产精品麻豆| 中文字幕高清在线视频| 综合色av麻豆| 国产激情偷乱视频一区二区| 十八禁网站免费在线| 精品国内亚洲2022精品成人| 亚洲av美国av| 美女被艹到高潮喷水动态| 国产精品一区二区三区四区久久| 可以在线观看毛片的网站| 亚洲av美国av| 国产精品电影一区二区三区| 久久久久久久午夜电影| 亚洲成av人片在线播放无| 18禁在线播放成人免费| 成人永久免费在线观看视频| 免费看美女性在线毛片视频| 757午夜福利合集在线观看| 久久久久久久精品吃奶| 黄色成人免费大全| 青草久久国产| 日韩免费av在线播放| 51国产日韩欧美| 97碰自拍视频| 又爽又黄无遮挡网站| 精华霜和精华液先用哪个| 国产乱人视频| 一级黄片播放器| 亚洲精品成人久久久久久| 国产av一区在线观看免费| 十八禁网站免费在线| 国产亚洲精品综合一区在线观看| 狠狠狠狠99中文字幕| 免费av不卡在线播放| 1024手机看黄色片| 无遮挡黄片免费观看| 亚洲久久久久久中文字幕| 88av欧美| 免费av观看视频| 狠狠狠狠99中文字幕| 久久香蕉精品热| 日韩有码中文字幕| 国产一区在线观看成人免费| 国产高清三级在线| 亚洲性夜色夜夜综合| 老司机在亚洲福利影院| 俺也久久电影网| 好男人电影高清在线观看| 中出人妻视频一区二区| 欧美黄色片欧美黄色片| 中文资源天堂在线| 一边摸一边抽搐一进一小说| 亚洲专区中文字幕在线| 欧美日韩亚洲国产一区二区在线观看| 五月玫瑰六月丁香| 日韩国内少妇激情av| 欧美+日韩+精品| 亚洲精品影视一区二区三区av| 精品一区二区三区视频在线观看免费| 日本撒尿小便嘘嘘汇集6| 色噜噜av男人的天堂激情| 亚洲在线自拍视频| 在线观看美女被高潮喷水网站 | av在线蜜桃| 老司机午夜十八禁免费视频| 99久久精品热视频| 波多野结衣巨乳人妻| 久久6这里有精品| 午夜福利在线观看免费完整高清在 | 欧美一级毛片孕妇| 美女黄网站色视频| 91九色精品人成在线观看| 日本 欧美在线| ponron亚洲| 亚洲精品456在线播放app | 一区福利在线观看| 免费人成在线观看视频色| 性色av乱码一区二区三区2| 久久香蕉精品热| 免费在线观看亚洲国产| 免费av毛片视频| 午夜免费男女啪啪视频观看 | 日本免费一区二区三区高清不卡| 99久久精品国产亚洲精品| 久久久成人免费电影| 丁香六月欧美| 久久亚洲真实| 久久精品人妻少妇| 成人18禁在线播放| 久久久久久大精品| 一级黄色大片毛片| 国产成人av教育| 亚洲美女视频黄频| 国产激情偷乱视频一区二区| 男女视频在线观看网站免费| 亚洲 国产 在线| 亚洲精品乱码久久久v下载方式 | 岛国在线观看网站| 午夜精品久久久久久毛片777| 精华霜和精华液先用哪个| 夜夜躁狠狠躁天天躁| 18禁美女被吸乳视频| 精品一区二区三区视频在线观看免费| 国产成人欧美在线观看| 国产高清视频在线观看网站| 国产精品乱码一区二三区的特点| 麻豆国产av国片精品| 成人三级黄色视频| 一级毛片女人18水好多| av视频在线观看入口| 国产高潮美女av| 国产中年淑女户外野战色| 此物有八面人人有两片| 桃色一区二区三区在线观看| 国产三级在线视频| 亚洲精品456在线播放app | 中文字幕久久专区| 国产一区二区激情短视频| 久久久久久久久中文| 国产精品女同一区二区软件 | 在线免费观看不下载黄p国产 | 免费人成在线观看视频色| 丰满乱子伦码专区| 国产免费男女视频| 亚洲黑人精品在线| 亚洲成人免费电影在线观看| 99在线视频只有这里精品首页| 又粗又爽又猛毛片免费看| 国产亚洲欧美98| 色噜噜av男人的天堂激情| 色尼玛亚洲综合影院| 久久香蕉国产精品| 变态另类成人亚洲欧美熟女| 夜夜夜夜夜久久久久| 五月玫瑰六月丁香| 99久久九九国产精品国产免费| 午夜精品久久久久久毛片777| av欧美777| a级一级毛片免费在线观看| 午夜福利在线观看吧| 日韩免费av在线播放| 国产精品永久免费网站| or卡值多少钱| 精品国产亚洲在线| 精品午夜福利视频在线观看一区| 成年女人永久免费观看视频| 小说图片视频综合网站| 国产单亲对白刺激| 有码 亚洲区| 麻豆成人av在线观看| 少妇人妻精品综合一区二区 | 叶爱在线成人免费视频播放| 国产极品精品免费视频能看的| 精品日产1卡2卡| 成人特级av手机在线观看| 中文亚洲av片在线观看爽| 中出人妻视频一区二区| 琪琪午夜伦伦电影理论片6080| 日本成人三级电影网站| 久久草成人影院| 国产国拍精品亚洲av在线观看 | 我要搜黄色片| 成年女人看的毛片在线观看| 欧美极品一区二区三区四区| or卡值多少钱| 不卡一级毛片| 亚洲国产精品久久男人天堂| 欧美一区二区亚洲| 免费在线观看影片大全网站| 欧美国产日韩亚洲一区| 观看美女的网站| 每晚都被弄得嗷嗷叫到高潮| 在线播放国产精品三级| 一个人免费在线观看电影| 精品电影一区二区在线| 天堂影院成人在线观看| 亚洲七黄色美女视频| 男人舔奶头视频| 淫妇啪啪啪对白视频| 嫩草影院入口| 亚洲欧美日韩无卡精品| 真人做人爱边吃奶动态| 青草久久国产| 搡女人真爽免费视频火全软件 | 淫秽高清视频在线观看| 国产99白浆流出| 久久伊人香网站| 日本免费a在线|