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

    Multi-objective steady-state optimization of two-chamber microbial fuel cells☆

    2017-05-29 10:47:56KeYangYijunHeZifengMa

    Ke Yang,Yijun He*,Zifeng Ma

    Shanghai Electrochemical Energy Devices Research Center,Department of Chemical Engineering,Shanghai Jiao Tong University,Shanghai 200240,China

    1.Introduction

    Microbial fuel cells(MFCs)provide a promising technology not only to renewable electricity generation,but also to wastewater treatment[1].MFCs can utilize a variety of materials as carbon sources such as natural organic matter,complex organic waste or renewable biomass.In comparison with the conventional fuel cells,MFCs use living organisms as a biocatalyst and often operate at mild conditions[2–5].Although MFCs have,in recent years,attracted significant attentions from both academia and industry,the intrinsic low biochemical reaction rate of MFCs technically hinders the successful commercial application of MFCs for energy generation or wastewater treatment.On the other hand,the complex bio-electrochemical reactions and transport phenomena are not well understood,which further limits the optimal design and operation ofMFCs.Hence,debottlenecking ofMFCs requires a breakthrough of multidisciplinary cooperation frommicrobiology,electrochemistry,electrical engineering and material science.

    Apart from experimental research,mathematical modeling approach provides an alternative promising tool to understand the reaction and transport mechanisms in MFCs.Both lumped and distributed parameter models with different complexities have been developed[6–13].Based on the developed models,itis able to investigate both the design parameters and operating condition effects on the performance of MFCs using sensitivity analysis,which helps to understand the importance ofparameters on system performance[11].Moreover,optimization method can be applied for attaining optimal design and operation of MFCs[13–16].Although the MFCs can simultaneously perform electricity generation and waste treatment,most of the published studies about performance optimization ofMFCs only focus on a single objective such as maximizing power output[13–16].Both experimental and simulation results have shown that the power output cannot be improved without the sacrifice of the maximal waste removal ratio.In addition,obtaining maximal power output often needs degrade the maximal attainable current density.That means these objectives,namely power density,attainable currentdensity and waste removalratio,mightbe conflicting,which implies that there does not exhibit a single solution that simultaneously optimizes each objective.Without a clear understanding of the relationships among those conflicting objectives,it is hard to make reasonable decisions for balancing trade-offs among conflicting objectives.In that case,multi-objective optimization methods are often used to obtain a set of Pareto-optimal solutions.Although it has been recognized that it is important to perform quantitative analysis of these conflicting objectives for obtaining a thorough understanding of their relationships,to our knowledge,there are rare studies about multi-objective optimization and multi-criteria decision making for optimal design and operation in the MFC community.In this study,we will try to establish a systematic multi-objective optimization framework to simultaneously maximize the power density,maximal attainable current density and waste removal ratio.The relationships among these objectives will be thoroughly investigated using the graphical visualization technique,i.e.level diagrams.In addition,the effects of operating conditions on those competing objective will be also investigated.

    The remainder of this paper is structured as follows.The mathematical model and multi-objective optimization problems of the twochamber MFCs are first introduced.The detailed implementations of the multi-objective genetic algorithm and the level diagrams method are next proposed.The results and discussion are then presented.Finally,conclusions and future research directions are provided.

    2.Two-chamber Microbial Fuel Cells

    2.1.Principle and experimental set-up of MFC

    MFCs can be broadly divided into two general categories,those that use a mediator and those that are mediator-less.These two types of MFCs are classified based on two transfer electron mechanisms,i.e.indirect electron transfer and direct electron transfer.The former uses a chemical mediator to transfer electron from microorganisms in the cell to the anodic electrode.The oxidized mediators arefirstly reduced by capturing electrons from within the membrane of and release the electrons to the electrode.Modern developed MFCs are often mediator-less,in which the microorganisms typically have electrochemically active redox proteins on their outer membrane that can transfer electrons directly to the anode.There are several electrochemically active microorganisms such asShewanella putrefaciens[17]andRhodoferax ferrireducens[18]which could form a bio film on the anode surface and transfer electrons directly to electrode across the membrane.These microorganisms often have high Coulombic efficiency and can avoid the use of the expensive and toxic mediators.

    A typical two-chamber MFC consists of an anodic chamber and a cathodic chamber separated by a cation-exchange membrane(CEM).In the anodic compartment,organic fuel is oxidized by microorganisms to generate carbon dioxide,electrons and protons.Electrons and protons are transferred to the cathodic compartment through an external electric circuit and the CEM,respectively,and then consumed at the cathode through combining with oxygen to form water.

    The schematic diagram of the experimental set-up of acetate two-chamber MFC can be referred to the studies reported elsewhere[13,16].The anode and cathode chambers are separated by a Na fion CEM.Acetate fuel and air-saturated water are continuously fed into the anode and cathode chambers respectively.Nitrogen gas is purged to keep the anode chamber anoxic.Each compartment contains two pieces of graphite felt as electrodes,but the graphite felt used for the cathode is coated with platinum powder.The multimeter is used to measure the potential between the anode and cathode,and a personal computer is used to record the potential through a data-acquisition system.Given a specific operating condition,the output power can be regulated through the change of the external resistance.

    2.2.Two-chamber MFC model

    Although there are different modeling approaches available in published literature concerning the MFCs,most reported models consisting of partial differential equations[7–11]have inherent complexity and heavy computational cost,which may not meetthe requirements of optimization and control in a real-time way and consequently may not be readily implemented by the majority ofthe MFC community.Recentdeveloped ordinary differential equations based lumped parameter models[12,13],which can be able to adequately capture MFC behavior atvarious operating conditions,could significantly decrease the computational cost and provide an alternative solution for meeting the requirements of optimal design,operation and control of MFCs.It is worth noting that although simplified models neglect spatial nonuniformities and cannot describe all the detailed phenomena occurring in the MFC,itmay be the only practicaland effective way for performing real-time optimization and control.Hence,in this study,the MFC model developed by Zenget al.[13]is used to perform multi-objective optimization.In this model,both the anodic and cathodic compartments are treated as a continuously stirred tank reactor(CSTR).It should be mentioned that although the CSTR model neglects the bio film formation on the anode side,both fitting and validation results have illustrated the effectiveness of this CSTR model,which is ready for further model analysis and optimization.The mathematical description of a two-chamber MFC CSTR model developed by Zenget al.[13]is briefly shown as follows.

    For acetate fueled two-chamber MFC,the oxidation of acetate and the reduction of dissolved oxygen are shown in Eqs.(1)and(2)at the anode and cathode,respectively.

    The rates of electrochemical reaction at the anode and cathode are formulated in Eqs.(3)and(4)using Butler–Volmer type expressions respectively.

    The mass balances of the four components in the anode chamber,namely acetate,dissolved carbon dioxide,hydrogen ion and biomass,are shown in Eqs.(5)–(8)respectively.

    where subscript‘a(chǎn)’denotes the anode.

    The mass balances of the three components in the cathode chamber,namely dissolved oxygen,hydroxyl and cation M+,are shown in Eqs.(9)–(11)respectively.

    where subscript ‘c’denotes the cathode;NMis the flux of M+ions transported from the anode to cathodeviathe CEM and is shown in Eq.(12).

    The charge balances at the anode and cathode are given by Eqs.(13)and(14)respectively.

    where ηaand ηcdenote the overpotential of anode and cathode respectively;Fis Faraday's constant andicellis the current density.

    The cellvoltageUcellis calculated as Eq.(15)by neglecting the ohmic drops in the current-collectors and electric connections.

    In practical,the output power is regulated through the change of external resistanceRex,and the cell current is calculated withIcell=Ucell/Rex.Note that the relationship between the current density and the cell current isicell=Icell/Am,whereAmis the area of membrane.Thus,the current density can be calculated as Eq.(16)by means of external resistance.

    By replacing the current density in Eqs.(12)–(14)with the external resistance using Eq.(16),the MFC model becomes more practical and could avoid determining the maximal current density for attaining meaningful system output.For solving the steady-state MFC model,ther1andr2can be firstly determined when using the current density as the input;the concentrations of seven state variables can be then calculated using Eqs.(5)–(11); finally,the overpotential of anode and cathode can be analytically computed using Eqs.(3)–(4)respectively.However,ifan inappropriate currentdensity is selected,the systemoutput would be illogical,e.g.the cell voltage mightbecome negative.If the external resistance is used,unreasonable system output can be effectively avoided.However,it is worth noting that the steady-state MFC model needs to be solved using a nonlinear optimization method because of the existence of coupling between external resistance and overpotential.

    The description and nominalvalues ofthe above parameters and operating conditions are summarized in Table 1 obtained from Zenget al.[13].Ithas been recognized thatthe performance ofMFCs is strongly dependent on operating conditions,such as pH,temperature,organic loading rate,feed rate and shear stress.Detailed analyses of the influences of operating conditions on the MFC performance can be referred to the reported review[2].In this study,four operating conditions,namely flow rate(Qa)and concentration(CiAnC)of acetate to anode,flow rate feeding to cathode(Qc)and external resistance(Rex)are selected as decision variables.

    2.3.Multi-objective optimization problem of MFCs

    A general mathematical formulation of multi-objective operational optimization problem for MFCs is as follows:

    where x and y denote the vectors ofcontrolvariables and state variables respectively;xminand xmaxare the lower and upper bounds of control variables respectively;Hrepresents the steady-state MFC modelobtained from Eqs.(3)–(16);F is theM-dimensionalobjective function vector,fm(x,y)indicates themth objective function.

    To evaluate the MFC performance,three objectives,namely power density,attainable current density and waste removal ratio,are used.

    Table 1Nominal values of parameters and operating conditions for MFC model[13]

    The power densityPDcellis calculated as

    Given specific values ofQa,CiAnCandQc,the currentdensity increases with the decreasing of the external resistance,and the maximal attainable current densityicell,maxcan be calculated using the minimal external resistanceRex,min.The larger value oficell,maxindicates the better current regulation ability of MFCs.

    The waste removal ratio WRR based on the inlet and outlet acetate concentration is defined as

    WRR provides a quantitative criterion of waste treatment ability.If the outlet acetate concentration is equal to 0,it implies that the overall waste is ideally removed and the value ofWRR is equalto 1;ifthe outlet acetate concentration is equal toCiAnC, it means that there is no waste is removed and the value of WRR is equal to 0.Note that the value of WRR lies in the range[0,1].The larger value of WRR indicates the better waste treatment ability of MFCs.

    Optimal operation of MFCs should balance the trade-offs among the following maximization criteria:power density,maximal attainable current density and waste removal ratio.It is worth noting that these objectives are generally conflicting,and consequently there exists no single solution that simultaneously optimizes all objective functions.Thatmeans none ofthe objectives can be improved withoutthe sacrifice of the other objectives.Therefore,it is desirable to use the multiobjective optimization method to generate a set of Pareto-optimal solutions for revealing the relationship among these conflict objectives.In this study,we will thoroughly investigate three bi-objective optimization problems,namely maximization of power density and waste removal ratio,maximization of power density and maximal attainable current density and maximization of waste removal ratio and maximal attainable current density,and one three-objective optimization problem,namely simultaneous maximization of power density,waste removal ratio and maximal attainable current density.

    3.Methodology

    3.1.Multi-objective optimization method

    Multi-objective optimization(MOO)is the process of simultaneously optimizing two or more conflicting objectives subject to certain constraints.MOO has been widely applied in many fields such as engineering and economics,in which optimal decisions need to be taken in the presence of trade-offs between two or more conflicting objectives.For a nontrivial MOO problem,there usually does not exhibit a single solution that simultaneously optimizes each objective and there exists a set of Pareto-optimal solutions.A solution is called nondominated or Pareto-optimal,if none of the objective functions can be improved in value without degrading some of the other objective values.The region defined by the objectives for all Pareto-optimal solutions is often called the Pareto front.The basic goal of MOO methods is to find a uniformly distributed set of Pareto-optimal solutions,which could approximate discretely the Pareto front.

    The MOO methods can be broadly decomposed into two categories:scalarization approaches and Pareto approaches.The former solves the MOOproblems by translating itback to a single(or a series ofsingle)objective scalarproblem.The formation ofaggregate objective function requires prescribinga prioripreferences or weights between objectives.Weighted sum approach,goal programming and lexicographic approach are three common types ofscalarization methods.In comparison with the scalarization approach,the Pareto approach could obtain a set of Pareto-optimal solutions in one run showing a more efficient way for solving MOO problems.Multi-objective genetic algorithms,adaptive weighted sum method and normal boundary intersection are three types of Pareto approach.In practice,it can be very difficult to precisely and accurately set up the preferences or weights,even for someone familiar with the problem domain.Therefore,the best way to solve the MOO problems is through the generation of a set of Pareto-optimal solutions using Pareto approaches,which could provide a wide range of decision options to the practitioners.

    Genetic algorithm(GA)is a kind of population-based optimization methods based on the principles of the evolutionviagenetics operation and natural selection and has been widely applied in engineering optimization[19].In a GA,a population of solutions to optimization problems evolves toward better solutions using techniques inspired by natural evolution,such as selection,mutation,and crossover.Non-dominated sorting genetic algorithm(NSGA)[20]and its enhanced variants,i.e.elitist NSGA-II[21]and controlled elitist NSGA-II[22],are an outstanding type of Pareto approach for solving MOO problems.The controlled elitist NSGA-II algorithm not only favors individuals with better fitness value,but also favors individuals that can help increase the diversity of the population even if they have a lower fitness value.The characteristics of trade-off between exploration and exploitation in controlled elitist NSGA-II could play an important role in maintaining the diversity of population for convergence to an optimal Pareto front.The details of controlled elitist NSGA-II algorithm can be referred to Deb and Goel's study[22]and the implementation procedure of NSGA-II is briefly described below.

    Step 1:Initially create a random parent populationP0of sizeNum,and assign each solution a fitness equal to its non-domination level obtained by fast non-dominated sorting method.

    Step 2:Use binary tournament selection,recombination and mutation operators to create a child populationQ0of sizeNum.

    Step 3:Sett=0.

    Step 4:Form a combined populationRt=Pt∪Qtof size 2Num,and each solution ofRtis assigned a fitness equal to its non-domination level.

    Step 5:Generate the new parentpopulationPt+1ofsizeNum,which uses a geometric distribution to maintain the number of individuals in each front.

    Step 6:Generate a new child population of sizeNumfrom parent populationPt+1by using binary tournament selection,recombination and mutation operators.

    Step 7:Sett=t+1;iftis less than maximum generation number,go back to Step 4;otherwise,exit and output the optimization results.

    3.2.Graphical visualization of Pareto front

    Once obtaining a set of Pareto-optimal solutions using NSGA-II,the next important step is to select one,or more,solutions inside the Pareto-optimal solutions.Without given sufficient preferences,visualization techniques could play a significant key role in helping decision-makers make proper decisions.Scatter diagrams and parallel coordinates are two common types of visualization techniques for multidimensional visualization of Pareto front.For the scatter diagrams,the complexity of the representation increases significantly with the dimension.The parallel coordinates,on the other hand,shows a compact representation through plotting a multidimensional point in a two-dimensional graph,but it would lose clarity and be difficult to perform analysis for large number of Pareto-optimal solutions.Level diagrams(LD)method,developed by Blascoet al.[23],is a new graphical visualization of both the Pareto front and set,and can be easily combined with preferences for helping the decisionmaker.The LD method firstly classifies the Pareto front points by means of their proximity to ideal points measured with a specific norm of normalized objectives and then synchronizes the objective and parameter diagrams.A description of the LD method is briefly described as follows.

    Let xi∈ X?? X denote theith Pareto-optimal solution,where X?is the obtained finite set of Pareto-optimal solutions using NSGA-II,andfm(xi) ∈ F(X?)denote themth objective value in Pareto front with respect to theith Pareto-optimal solution.The maximum and minimum values of themth objective on the approximated Pareto front are given by Eqs.(20)and(21),respectively.

    A smaller value of‖F(xiàn)(x)‖p(p=1,2,∞)indicates a better solution.Based on the above normalization and norm computation procedures,every Pareto-optimal solution is classified and the points of the Pareto front are sorted in ascending order of the value of‖F(xiàn)(x)‖p.In the level diagrams,each objective or decision variable has its own graphical representation.The verticalaxis on allgraphs corresponds the levelvalue of‖F(xiàn)(x)‖p,which means that all graphs are synchronized with respect to this axis.The horizontal axis corresponds to values of the objectives,or decision variables.Note that level diagrams draw a point at the same level for all graphs in objective and decision variable space,which is useful for making proper decisions.Further detailed information about the LD method can be referred to the study of Blascoet al.[23].It should be noted that although the classical scatter diagrams are sufficient for a bi-objective problem,the level diagrams can provide more information about the distance to the ideal point,which is useful for decision-maker.In this study,the LD method is only carried out for the three-objective optimization problem.

    4.Results and Discussion

    The above described controlled elitist NSGA-II and level diagrams are applied to multi-objective operational optimization of twochamber MFC.Three bi-objective optimization and one threeobjective optimization problem are investigated.Both controlled elitist NSGA-II algorithm and LD method are implemented in MATLAB,and all computations are carried out on a PC with 2.60 GHz processor and 8 GB of RAM.The parameter values of controlled elitist NSGA-II algorithm are shown as follows:number of generations is 500,population size is 500,cross fraction is 0.8,uniform mutation rate is 0.01,and Pareto fraction is 0.25.The lower bounds ofQa,CiAnC,QcandRexare 0.5×10-5m3·h-1,1 mol·m-3,0.5 × 10-3m3·h-1and 10 Ω,respectively;the upper bounds ofQa,CiAnC,QcandRexare 3 × 10-5m3·h-1,2 mol·m-3,2 × 10-3m3·h-1and 10000 Ω,respectively.

    4.1.Maximizing power density and waste removal ratio

    The Pareto front obtained by maximizing power density and waste removal ratio without considering the criterion of attainable current density is presented in Fig.1.The Pareto front is smooth and well distributed over a wide range.It is found that as the power density increases,the waste removal ratio would decrease,indicating that a conflict exists between power density and waste removal ratio.It also means that improvement of power density would result in deterioration of waste removal ratio,and vice versa.From Fig.1,the maximum power density of 3.82 W·m-2corresponds to the minimum waste removalratio of0.67,while the maximumwaste removalratio of0.97 corresponds to the minimum power density of 0.64 W·m-2.It can be observed that at high power density region,waste removal ratio can be significantly improved with a relatively small sacrifice in power density.For example,increasing of waste removal ratio from 0.67 to 0.80 only leads to decreasing of power density from 3.82 to 3.50 W·m-2.While at high waste removal ratio region,achieving a relatively small improvement of waste removal ratio needs a much higher sacrifice in power density.For example,increasing of waste removal ratio from 0.90 to 0.97 would result in a significant decreasing of power density from 2.52 to 0.64 W·m-2.If a MFC is designed to only generate electricity,the Pareto point with maximum power density of 3.82 W·m-2becomes a best option.On the other hand,if a MFC is designed to wastewater treatment,the Pareto point with maximum waste removal ratio of 0.97 is a best option.However,practical MFC operation often avoids the extreme objective values and calls for a trade-off between power density and waste removalratio.Moreover,although the obtained Pareto front provides a wide range of competing options for optimal operation from a multi-objective perspective,additional preferences should be provided for assisting in final decision-making.The maximal attainable currentdensity(icell,max)againstpowerdensity and waste removal ratio corresponding to the Pareto front in Fig.1 is also plotted in Fig.2.Itis observed that(1)the maximum oficell,maxis 8.87 A·m-2;and(2)the maximal attainable current density increases with the increasing of power density,but with the decreasing of waste removal ratio.It is implied that without considering the criterion of maximization oficell,max,maximal attainable current density and power density seem to be not conflicting,where improving power density would also result in increasing attainable current density,whereas maximal attainable currentdensity and waste removalratio are conflicting,where waste removal ratio is improved with the sacrifice of attainable current density.However,itshould be noted thatiftreating attainable currentdensity as optimization objective,the exhibition of a clear conflict between power density and attainable current density would be illustrated and the maximum oficell,maxcan be significantly improved.

    Fig.1.Pareto front obtained from the simultaneous maximization of power density and waste removal ratio.

    To quantify the decision variables'role in realizing the competing objective values shown in Fig.1,Fig.3 plots the decision variables against power density and waste removal ratio.Fig.3(A)shows a clear increase in flow rate of fuel feed to anode(Qa)with increasing power density and decreasing waste removal ratio.The possible reason is that increasing ofQaresults in an increase of available carbon source,which favors to improve power output.However,it should be noted that the maximum value ofQais 1.1 × 10-5m3·h-1corresponding to maximal power density of 3.82 W·m-2in Fig.1,which only lies in the medium region in the search space.It is,thus,implied that further increasing ofQawill not improve power output.On the other hand,increasing ofQawould decrease the residence time of reactants,which might result in reducing waste removal capability.Fig.3(B)shows that almost the largest available concentration of fuel(CiAnC) is chosen uniformly for the entire Pareto.It indicates that increasing concentration of fuel is greatly helpful to increase both power density and waste removal ratio.Fig.3(C)indicates that the values of flow rate of feeding to cathode(Qc)are scattering and have no clear trends compared to the smooth variation of the Pareto front they generated(Fig.1).A two-chamber MFC is often operated with an excessive oxygen supply in the cathode for maximizing power output and stimulating waste removal.Based on our simulation results,it is found that if fixing a relatively highQcvalue of 2.0×10-3,the obtained Pareto front is almost the same as thatin Fig.1.Itis,thus,implied that maintaining a technical feasible high enough flow rate of feeding to cathode(Qc)could be helpful to maximize power output and waste removal.However,it is worth noting that if a MFC has been operated at an anode reaction limiting region,a furtherincrease ofQcwould be meaninglessand cannotimprove power output and waste removal.From Fig.3(D),it can be observed that(1)at the low power density region of[0.60,1.44]W·m-2and high waste removal ratio region of[0.95,0.97],external resistance(Rex)approximately increases linearly from 70 to 170 Ω with the increasing of power density and the decreasing of waste removal rate;and(2)at the power density region of[1.44,3.82]W·m-2and waste removal ratio region of[0.64,0.95],the values of internal resistance(Rex)almost uniformly scatter at the region of[170,190]Ω.In general,the external resistance is optimally manipulated to regulate the current for tracking the maximal power point.Fig.4 plots the currentdensity(icell)againstexternalresistance corresponding to the Pareto frontin Fig.1.When the externalresistance ismanipulated from70 to 170 Ω,the current density of MFC almost remains constant at 4.2 A·m-2.When the external resistance varies from 170 to 190 Ω,the values of current density scatter at the region of 4.2-6.6 A·m-2.Moreover,the currentdensity againstpowerdensity and waste removal ratio is also plotted in Fig.5.Itis observed that increasing of power density mostly requires the increasing of current density,whereas increasing of waste removal ratio is mostly attributed to the decreasing of current density.

    Fig.2.Maximal attainable current density versus power density()and waste removal ratio(○)corresponding to the Pareto front in Fig.1.

    Fig.3.Decision variables versus power density()and waste removal ratio(○)corresponding to the Pareto front in Fig.1.

    Fig.4.Externalresistance versus currentdensity corresponding to the Pareto front in Fig.1.

    Fig.5.Current density versus power density()and waste removal ratio(○)corresponding to the Pareto front in Fig.1.

    4.2.Maximizing power density and maximal attainable current density

    Fig.6 shows the Pareto frontby maximizing power density and attainable current density without considering the criterion of waste removal ratio.It is observed that power density decreases with the increasing of the maximal attainable current density.The maximum power density of 3.78 W·m-2corresponds to the minimum maximal attainable current density of 9.96 A·m-2,while the maximum maximal attainable current density of 14.88 A·m-2corresponds to the minimum power density of 2.28 W·m-2.In comparison with Fig.2,it is obviously found that maximalattainable currentdensity can be significantly improved when explicitly considering it as an optimization objective.Moreover,it is observed that at high maximal attainable current density,achieving a relatively small improvement of maximal attainable current density needs a much higher sacrifice in power density.For example,increasing of maximal attainable currentdensity from14.49 to 14.88 A·m-2would resultin a significant decreasing of power density from 2.80 to 2.28 W·m-2.Fig.7 plots the waste removal ratio against power density and maximal attainable current density corresponding to the Pareto front in Fig.6.It is found that waste removal ratio increases with the increasing of power density,whereas it increases with the decreasing of maximal attainable current density.The maximum waste removal ratio is only 0.67,which approximately corresponds to the minimum waste removal ratio in Fig.1.This means withoutexplicitly considering waste removalratio as optimization objective,itfailed to obtain a relatively high waste removalratio.Hence,it is important to explicitly introduce waste removal ratio as an optimization objective if a MFC is designed to simultaneous electricity generation and waste treatment.

    Fig.8 plots the decision variables against power density and maximal attainable current density corresponding to the Pareto front in Fig.6.From Fig.8(A),it is observed that power density decreases with the increasing ofQa,whereas maximal attainable current density increases with the increasing ofQa.It is also found that the operating region ofQais(1.3-2.6)× 10-5m3·h-1,which is significantly higher than that in Fig.3(A).Fig.8(B)shows that almost the largest availableCiAnCof 2 mol·m-3is chosen for the entire Pareto.It is thus indicated that achieving higher maximal attainable current density requires a substantial increase ofQa,but increasing ofQacannot further improve the maximal attainable current density without increasingCiAnC. Fig.8(C)shows that almostQcscatters at a narrow region of(1.95-2.0)× 10-3m3·h-1for the entire Pareto.Fig.8(D)shows that at the low power density region of 2.28-3.20 W·m-2,almostRexis chosen in a narrow region of 76-82 Ω,whereas at the high power density region of 3.20-3.78 W·m-2,power density approximately linearly increases with the increasing ofRex.Note that the plot of maximal attainable current density againstRexis meaningless,because maximal attainable current density is calculated at the minimal value ofRexof 10 Ω.Fig.9 plots the current density(icell)againstRexcorresponding to the Pareto front in Fig.6.WhenRexis manipulated from 76 to 90 Ω,the values oficellscatter at the region of 7.7-9.2 A·m-2.WhenRexvaries d from 90 to 142Ω,icellshows an approximately lineardecreasing trend with the increasing ofRex.Moreover,the current density against power density and maximal attainable current density is plotted in Fig.10.It is observed that at the high power density region of 3.20-3.78 W·m-2,power density decreases with the increasing oficell,whereas at the low power density region of 2.28-3.20 W·m-2,power density increases with the increasing oficell.

    Fig.6.Pareto front obtained from the simultaneous maximization of power density and maximal attainable current density.

    Fig.7.Waste removal ratio versus power density()and maximal attainable current density(○)corresponding to the Pareto front in Fig.6.

    4.3.Maximizing waste removal ratio and maximal attainable current density

    Fig.11 shows the Pareto front by maximizing waste removal ratio and attainable current density without considering the criterion of power density.It is observed that waste removal ratio decreases with the increasing of the maximalattainable currentdensity.The maximum waste removal ratio of 0.98 corresponds to the minimum maximal attainable current density of 4.20 A·m-2,while the maximum maximal attainable currentdensity of14.88 A·m-2corresponds to the minimum waste removalratio of0.67.Atthe high maximalattainable currentdensity region of 12.96-14.88 A·m-2,a relatively small improvement oficell,maxwould result in a significant deterioration of waste removal ratio.Fig.12 plots the power density against waste removal ratio and maximal attainable current density corresponding to the Pareto front in Fig.11.It is observed that power density approximately linearly increases with the increasing of maximal attainable current density,whereas it decreases with the increasing of waste removal ratio.It is also found that the maximum power density is only 1.16 W·m-2without explicitly considering it as an optimization objective.It is worth noting that the indication of maximal attainable current density is considered to be of importance only when it is combined with the power density for evaluating the MFC performance in electrical applications.This study only attempts to illustrate the relationship between maximal attainable current density and waste removal ratio,practical optimization of maximal attainable current density should be performed simultaneously with power density.

    Fig.9.Externalresistance versus currentdensity corresponding to the Pareto front in Fig.6.

    Fig.10.Current density versus power density()and maximal attainable current density(○)corresponding to the Pareto front in Fig.6.

    Fig.11.Pareto frontobtained from the simultaneous maximization of waste removalratio and maximal attainable current density.

    Fig.12.Power density versus waste removal ratio()and maximal attainable current density(○)corresponding to the Pareto front in Fig.11.

    Fig.13 plots the decision variables against waste removal ratio and maximal attainable current density corresponding to the Pareto front in Fig.11.It is observed from Fig.13(A)that waste removal ratio decreases with the increasing ofQa,whereas maximal attainable current density increases with the increasing ofQa.The operating region ofQais(0.5-2.6) × 10-5m3·h-1,which is much wider than that both in Figs.3(A)and 8(A).From Fig.13(B),it is observed that almost the largest availableCiAnCof 2 mol·m-3is chosen for the entire Pareto,which is similar with Figs.3(B)and 8(B).Fig.13(C)and(D)shows that almost the values ofQcscatter at the region of 1.3-2.0 m-3·h-1and almost the values ofRexscatter at a narrow region of 11.1-11.5 Ω.From Fig.14,it is found that the values oficellcan vary from 4.2 to 14.9 A·m-2with a small variation inRex.The large variation inicellis mainly attributed to the large variation inQa.Fig.15 shows the plot oficellagainst maximal attainable current density and waste removal ratio corresponding to the Pareto front in Fig.11.It is observed that waste removal ratio shows a decreasing trend with the increasing oficell.The linear relationship with an approximate slope of 1 between maximal attainable current density andicellis attributed to the fact that almost the values of Pareto optimalRexapproach its minimum value of 10 Ω,at which the maximal attainable current density is calculated.Moreover,from Figs.12 and 15,it is found that power density approximately increases with the increasing oficell.

    4.4.Maximizing power density,waste removal ratio and maximal attainable current density

    Fig.13.Decision variables versus waste removal ratio()and maximal attainable current density(○)corresponding to the Pareto front in Fig.11.

    The Pareto front obtained by simultaneously maximizing power density,waste removal ratio and maximal attainable current density is shown in Fig.16.It is obtained that(1)the maximum power density of 3.82 W·m-2corresponds to the maximal attainable current density of 8.86 A·m-2and the waste removal ratio of 0.67;(2)the maximum waste removal ratio of 0.97 corresponds to the power density of 0.39 W·m-2and the maximal attainable current density of 4.21 A·m-2;and(3)the maximum maximal attainable current density of 14.88 A·m-2corresponds to the power density of 2.28 W·m-2and the waste removalratio of0.35.Moreover,the minima ofpower density,waste removal ratio and maximal attainable current density are 0.39 W·m-2,0.35 and 4.10 A·m-2,respectively.It is thus implied that achieving the highest waste removal ratio would result in the minimum power output,whereas a highest maximal attainable current density is achieved with the lowest waste removal ratio.

    Fig.17 shows a 2-Dvisualization ofthe Pareto frontby the projection ofthe 3-DPareto frontin Fig.16.Itis observed from Fig.17(A)that(1)at the high waste removal ratio region of 0.80-0.97,it decreases with the increasing of power density;(2)at the low waste removal region of 0.35-0.55,it increases with the increasing of power density;and(3)at the medium waste removal region of 0.55-0.80,the relationship between waste removal ratio and power density shows scattering characteristic.Moreover,it is observed from Fig.17(A)that the low power density region of 0.39-2.30 W·m-2is mostly obtained at the high waste removal ratio region and the high power density region of 3.0-3.82 W·m-2is obtained at the medium waste removal ratio region.From Fig.17(B),it can be seen that(1)at the low maximal attainable current density region of 4.10-9 A·m-2,power density approximately increases with the increasing of maximal attainable current density;(2)at the high maximal attainable current density region of 12-15 A·m-2,increasing of power density would result in a decreasing of maximal attainable current density;and(3)at the medium maximal attainable current density region of 9-12 A·m-2,the changes in power density with respect to maximal attainable current density show scattering characteristics.It is seen from Fig.17(C)that at both the low and high waste removal ratio regions,waste removal ratio decreases with the increasing of maximal attainable current density,whereas at the medium waste removal ratio region,waste removal ratio shows a decreasing scattering characteristic with respect to maximal attainable current density.

    According to the Pareto front in Fig.16,the LD method is applied to aid in graphical visualization and optimal decision making.Fig.18 shows the 2-norm level diagram representation of the Pareto front in Fig.16.It is found that the lowest level value is always over 0.52.That means the Pareto front is relatively far from the ideal point.It should be noted that in the LD method,all points have been normalized and the ideal point is thus of[0,0,0]for three-objective maximization problems.Moreover,note that the points situated at the lower levels correspond to the zones of the Pareto front nearer to the ideal point.It is observed from Fig.18 that there exhibits several points at the lower levels,which means that there are several points relatively close to the ideal point.The nearest points to the ideal are approximately at 2.84-3.25 W·m-2for power density,0.69-0.78 for waste removal ratio and 10.82-12.19 A·m-2for maximal attainable current density.These points provide a set of promising candidates for final optimal decision making.However,it should be noted that the LD method with 2-norm treats each objective equally,and level computation is based on the distance from the ideal point.The final decision would be significantly affected if each criterion is given a proper preference or constraint.If power density needs to be above 3.50 W·m-2,it can be seen from Fig.18 that the 2-norm level should be at the range of 0.57-0.74.At this level,the waste removal ratio and maximal attainable current density are situated at the ranges of 0.57-0.73 and 8.73-11.81 A·m-2,respectively.If waste removal ratio needs to be above 0.95,it is found that the 2-norm level should be above 1.20,which means that these points are far from the ideal point.At this level,the power density and maximal attainable current density are situated at the ranges of 0.39-1.09 W·m-2and 4.10-5.88 A·m-2,respectively.That means achieving extreme high waste removal ratio needs significantly sacrifice power density and maximal attainable current density.Moreover,if maximal attainable current density needs to be above 14.0 A·m-2,the 2-norm level should be at the range of 0.74-1.10,which means that these points are relatively far from the ideal point.At this level,the power density and waste removal ratio are situated at the ranges of 2.28-3.03 W·m-2and 0.35-0.55,respectively.It also implies that achieving the highest maximal attainable current density would result in the lowest waste removal ratio.

    Fig.14.External resistance versus current density corresponding to the Pareto front in Fig.11.

    Fig.15.Current density versus waste removal ratio()and maximal attainable current density(○)corresponding to the Pareto front in Fig.11.

    Fig.16.Pareto front obtained from the simultaneous maximization of power density,waste removal ratio and maximal attainable current density.

    Fig.19 shows the 2-norm level diagrams representation of the decision variables corresponding to the Pareto front in Fig.16.It is observed that(1)the values ofQaare situated at the range of(0.5-2.56)× 10-5m3·h-1;(2)the values ofQcare almost over 1.5× 10-5m3·h-1;(3)the values ofCiAnCmostly reach its highest limit of 2 mol·m-3,which is in accordance with the three bi-objective optimization results;(4)the values ofRexare below 200 Ω and most ofRexare manipulated at the range of 50-80Ω,and(5)the values oficellare regulated atthe range of 4.07-10.56 A·m-2.Moreover,it is found that the ranges of operating variables corresponding to the lower level are approximately atQa∈ [1.4,1.7]× 10-5m3·h-1,Qc∈ [1.5,1.8]× 10-3m3·h-1,CiAnC≈ 2 mol·m-3,Rex∈ [50,70]Ω andicell∈[9.44,10.56]A·m-2.

    Note that Blascoet al.have recommended representing the same Pareto front with different norms to see the difference[23].However,it is found that for this three objective optimization problem,both 1-norm and in finity norm representations do not illustrate much new information by comparing to 2-norm representation,and we neglect to depict these two representations.Moreover,although the level diagrams provide a promising graphical analysis tool for high-dimensional Pareto front,more sophisticated preferences need to be introduced for final decision making.

    Fig.17.2-D visualization of Pareto front by the projection of a 3-D Pareto front in Fig.16:(A)power density versus waste removal ratio;(B)power density versus maximal attainable current density;(C)waste removal ratio versus maximal attainable current density.

    5.Conclusions

    Multi-objective steady-state optimization of MFC was successfully carried out to understand the complex relationship between power density,waste removal ration and maximal attainable current density.Three bi-objective and one three-objective optimization problems were thoroughly investigated and the LD method was applied to assist in graphical analysis and decision making for three-objective optimization problem.The main findings of the present study can be summarized as follows:

    Fig.18.2-norm level diagram representation of the Pareto front in Fig.16.

    (1)Results ofthree bi-objective optimization problems indicated that power density,waste removal ratio and maximal attainable current density are mutually conflicting.The maximum values of power density,waste removal ratio and maximal attainable current density are 3.82 W·m-2,0.97 and 14.88 A·m-2,respectively.The effects offour operating conditions,namelyQa,Qc,CiAnCandRex,on MFC performance illustrated that(a)the flow rate of feeding to cathode(Qc)can be usually chosen to supply excessive oxygen;(b)the flow rate of fuel feed to anode(Qa)plays significant and different effects on objectives in different bi-objective problems and should be properly determined for balancing tradeoffs between conflicting objectives;(c)the acetate concentration(CiAnC)often reaches its highest limit of 2 mol·m-3for all problems;and(d)the externalresistance(Rex)should be properly manipulated below 200 Ω for different problems.

    (2)The three-objective optimization results demonstrated that these three objectives are conflicting in nature.Graphicalanalysis ofPareto front using level diagrams illustrated that the nearest points to the ideal point are at 2.84-3.25 W·m-2for power density,0.69-0.78 for waste removal ratio and 10.82-12.19 A·m-2for maximal attainable current density.Moreover,it was found that(a)achieving the high power density range of 3.50-3.82 W·m-2would result in the medium waste removal ratio range of 0.57-0.73 and maximal attainable current density range of 8.73-11.81 A·m-2;(b)achieving the high waste removal ratio above 0.95 would lead to low power density below 1.09 W·m-2and maximal attainable current density below 5.88 A·m-2;and(c)achieving the high maximal attainable current density above 14.0 A·m-2would resultin low waste removal ratio below 0.55.

    Hence,the aforementioned findings indicate that the integrated methodology of multi-objective genetic algorithm and level diagram providesa promising approach to trade-offconflicting objectives foroptimal operation of MFC.However,the present study only investigates the steady-state operational optimization of MFC.For improving the dynamic operational performance,future work should focus on multiobjective dynamic optimization of MFC.

    [1]B.E.Logan,Microbial Fuel Cells,Wiley,Hooken,NJ,2008.

    [2]V.B.Oliveirab,M.Sim?es,L.F.Melo,A.M.F.R.Pinto,Overview on the developments of microbial fuel cells,Biochem.Eng.J.73(2013)53–64.

    [3]B.E.Logan,B.Hamelers,R.Rozendal,U.Schr?der,J.Keller,S.Freguia,P.Aelterman,W.Verstraete,K.Rabaey,Microbial fuel cells:Methodology and technology,Environ.Sci.Technol.40(17)(2006)5181–5192.

    [4]B.H.Kim,H.J.KIM,M.S.Hyun,D.Park,Direct electrode reaction of Fe(III)-reducing bacterium,Shewanella putrefaciens,J.Microbiol.Biotechnol.9(2)(1999)127–131.

    [5]E.Katz,A.N.Shipway,I.Wilner,in:W.Vielstich,A.Lamm,H.A.Gasteiger(Eds.),Handbook of Fuel Cells:Fundamentals,Technology,Application,Wiley,Chichester,United Kingdom,2003.

    [6]X.C.Zhang,A.Halme,Modelling of a microbial fuel cell process,Biotechnol.Lett.17(2)(1995)809–814.

    [7]A.K.Marcus,C.I.Torres,B.E.Rittmann,Conduction-based modeling of the bio film anode of a microbial fuel cell,Biotechnol.Bioeng.98(2007)1171–1182.

    [8]C.Picioreanu,I.M.Head,K.P.Katuri,M.C.M.van Loosdrecht,K.Scott,Acomputational model for bio film-based microbial fuel cells,Water Res.41(3)(2007)2921–2940.

    [9]C.Picioreanu,I.M.Head,K.P.Katuri,M.C.M.van Loosdrecht,K.Scott,Mathematical model for microbial fuel cells with anodic bio films and anaerobic digestion,Water Sci.Technol.57(2008)965–971.

    [10]C.Picioreanu,K.P.Katuri,M.C.M.van Loosdrecht,I.M.Head,K.Scott,Modelling microbial fuel cells with suspended cells and added electron transfer mediator,J.Appl.Electrochem.40(2010)151–162.

    [11]C.Picioreanu,M.C.M.van Loosdrecht,T.P.Curtis,K.Scott,Model based evaluation of the effect of pH and electrode geometry on microbial fuel cell performance,Bioelectrochemistry78(1)(2010)8–24.

    [12]R.P.Pinto,B.Srinivasan,M.-F.Manuel,B.Tartakovsky,A two-population bioelectrochemical model of a microbial fuel cell,Bioresour.Technol.101(14)(2010)5256–5265.

    [13]Y.Z.Zeng,Y.F.Choo,B.H.Kim,P.Wu,Modelling and simulation of two-chamber microbial fuel cell,J.Power Sources195(1)(2010)79–89.

    [14]L.Woodward,M.Perrier,B.Srinivasan,R.P.Pinto,B.Tartakovsky,Comparison of real-time methods for maximizing power output in microbial fuel cells,AIChE J.56(10)(2010)2742–2750.

    [15]R.P.Pinto,B.Srinivasan,S.R.Guiot,B.Tartakovsky,The effect of real-time external resistance optimization on microbial fuel cell performance,Water Res.45(4)(2011)1571–1578.

    [16]Y.J.He,Z.F.Ma,Robust optimal operation of two-chamber microbial fuel cell system under uncertainty:a stochastic simulation based multi-objective genetic algorithm approach,Fuel Cells13(3)(2013)321–335.

    [17]H.J.Kim,H.S.Park,M.S.Hyun,I.S.Chang,M.Kim,B.H.Kim,A mediator-less microbial fuel cell using a metal reducing bacterium,Shewanella putrefaciens,Enzyme Microb.Technol.30(2)(2002)145–152.

    [18]S.K.Chaudhuri,D.R.Lovley,Electricity generation by direct oxidation of glucose in mediatorless microbial fuel cells,Nat.Biotechnol.21(2003)1229–1232.

    [19]M.Gen,R.Cheng,Genetic Algorithms and Engineering Optimization,Wiley-Interscience,Hooken,NJ,USA,1999.

    [20]N.Srinivas,K.Deb,Muiltiobjective optimization using nondominated sorting in genetic algorithms,Evol.Comput.2(3)(1994)221–248.

    [21]K.Deb,S.Agrawal,A.Pratap,T.Meyarivan,A fast elitist non-dominated sorting genetic algorithm for multiobjective optimization:NSGA-II,Proceedings of the Parallel Problem Solving From Nature VI Conference,Paris,France,2000.

    [22]K.Deb,T.Goel,Controlled elitist non-dominated sorting genetic algorithm for better convergence,Proceeding EMO'01 Proceedings of thefirst International Conference on Evolutionary Multi-criterion Optimization,Springer-Verlag,London,United Kingdom,2001.

    [23]X.Blasco,J.M.Herrero,J.Sanchis,M.Martínez,A new graphical visualization of ndimensional Pareto front for decision-making in multiobjective optimization,Inf.Sci.178(20)(2008)3908–3924.

    午夜福利,免费看| 六月丁香七月| 夜夜看夜夜爽夜夜摸| 一级毛片黄色毛片免费观看视频| 亚洲电影在线观看av| 国产精品一区二区在线观看99| 夫妻性生交免费视频一级片| 久久鲁丝午夜福利片| 免费看光身美女| 女人久久www免费人成看片| h日本视频在线播放| 国产精品国产三级专区第一集| 久久久久久久久大av| av国产精品久久久久影院| 97精品久久久久久久久久精品| 国产精品麻豆人妻色哟哟久久| 成人亚洲精品一区在线观看| 精品国产一区二区久久| av国产精品久久久久影院| 日韩一区二区三区影片| 黄色欧美视频在线观看| 国产精品.久久久| 日韩av不卡免费在线播放| 精品久久久久久电影网| 成人无遮挡网站| 久久久久久久大尺度免费视频| 国产毛片在线视频| 最黄视频免费看| 成人影院久久| 国产精品国产三级国产专区5o| 秋霞伦理黄片| 色网站视频免费| 久久精品国产亚洲av涩爱| 亚洲av日韩在线播放| 婷婷色综合大香蕉| 久久久久久久精品精品| 国产精品国产三级国产专区5o| 国产在线男女| 一个人免费看片子| 性高湖久久久久久久久免费观看| 黄片无遮挡物在线观看| 在线观看国产h片| 国产 精品1| 日本欧美视频一区| 中文天堂在线官网| 午夜激情久久久久久久| 乱系列少妇在线播放| 久久精品国产a三级三级三级| 亚洲久久久国产精品| 男女边摸边吃奶| 亚洲国产欧美在线一区| 一本久久精品| 熟妇人妻不卡中文字幕| 一级毛片黄色毛片免费观看视频| 国产成人免费观看mmmm| 精品久久国产蜜桃| 妹子高潮喷水视频| 国产精品.久久久| 婷婷色麻豆天堂久久| 多毛熟女@视频| 久久久久精品性色| 日日摸夜夜添夜夜添av毛片| 人妻夜夜爽99麻豆av| 日韩中文字幕视频在线看片| 熟女人妻精品中文字幕| 人妻一区二区av| 欧美亚洲 丝袜 人妻 在线| 伦理电影免费视频| 亚洲国产欧美日韩在线播放 | 在线观看美女被高潮喷水网站| 国产精品一区www在线观看| 亚洲av在线观看美女高潮| 麻豆成人av视频| 国产成人精品久久久久久| 黄片无遮挡物在线观看| 一级二级三级毛片免费看| 在现免费观看毛片| 免费人成在线观看视频色| 国产精品一区二区性色av| 自拍欧美九色日韩亚洲蝌蚪91 | 国产极品天堂在线| 欧美精品国产亚洲| 国产视频内射| av有码第一页| 黑人高潮一二区| 国产精品久久久久久av不卡| 成年av动漫网址| 久久久久国产网址| 亚洲国产av新网站| 色婷婷久久久亚洲欧美| 国产精品人妻久久久久久| 最新的欧美精品一区二区| 青春草亚洲视频在线观看| 免费在线观看成人毛片| 亚洲中文av在线| 免费播放大片免费观看视频在线观看| 国产亚洲精品久久久com| 青春草视频在线免费观看| 亚洲欧美精品自产自拍| 多毛熟女@视频| 狂野欧美白嫩少妇大欣赏| 国产精品人妻久久久久久| 免费高清在线观看视频在线观看| 少妇猛男粗大的猛烈进出视频| 国产无遮挡羞羞视频在线观看| 人人妻人人看人人澡| 久久久久久久久久久久大奶| freevideosex欧美| 国产伦精品一区二区三区视频9| 青春草亚洲视频在线观看| 亚洲av免费高清在线观看| 成人综合一区亚洲| 国产精品一区二区性色av| 免费av中文字幕在线| 亚洲四区av| 天堂俺去俺来也www色官网| 午夜福利,免费看| 亚洲精品自拍成人| 国产亚洲最大av| 日本猛色少妇xxxxx猛交久久| 91久久精品国产一区二区成人| 欧美亚洲 丝袜 人妻 在线| 久久婷婷青草| 国产伦精品一区二区三区四那| 亚洲,欧美,日韩| 精品国产一区二区三区久久久樱花| 寂寞人妻少妇视频99o| 久热这里只有精品99| 另类亚洲欧美激情| 久久精品国产亚洲av涩爱| 美女脱内裤让男人舔精品视频| 成年人免费黄色播放视频 | 在线观看免费日韩欧美大片 | 亚洲精品日韩av片在线观看| 日韩中字成人| 亚洲欧美成人精品一区二区| 国产精品久久久久久精品电影小说| av不卡在线播放| 国产精品久久久久久精品电影小说| 秋霞在线观看毛片| 中文字幕制服av| 精品人妻一区二区三区麻豆| 亚洲第一av免费看| a级片在线免费高清观看视频| 黄色日韩在线| 免费观看av网站的网址| 免费观看av网站的网址| 精华霜和精华液先用哪个| 精品一区在线观看国产| 在线 av 中文字幕| 80岁老熟妇乱子伦牲交| 99热网站在线观看| av在线播放精品| 亚洲成人一二三区av| 伦精品一区二区三区| 中文欧美无线码| 一级毛片电影观看| 美女内射精品一级片tv| 国产爽快片一区二区三区| 欧美日本中文国产一区发布| 国产国拍精品亚洲av在线观看| 在现免费观看毛片| 日韩,欧美,国产一区二区三区| 交换朋友夫妻互换小说| av.在线天堂| 永久网站在线| 欧美激情极品国产一区二区三区 | av网站免费在线观看视频| 国产极品天堂在线| 男女啪啪激烈高潮av片| 亚洲av国产av综合av卡| 男女边摸边吃奶| freevideosex欧美| 欧美精品国产亚洲| 美女xxoo啪啪120秒动态图| 最近中文字幕高清免费大全6| 一区二区三区四区激情视频| 久久久久久久精品精品| 成人二区视频| 又爽又黄a免费视频| 大陆偷拍与自拍| 日本猛色少妇xxxxx猛交久久| 亚洲欧美一区二区三区国产| 精品视频人人做人人爽| 国内揄拍国产精品人妻在线| 美女cb高潮喷水在线观看| 在线观看免费视频网站a站| 亚洲av中文av极速乱| 日韩成人av中文字幕在线观看| 韩国高清视频一区二区三区| 成人免费观看视频高清| 女人精品久久久久毛片| 欧美精品亚洲一区二区| 国产在线免费精品| 我的女老师完整版在线观看| 丝袜在线中文字幕| 熟女电影av网| 校园人妻丝袜中文字幕| av黄色大香蕉| 精品熟女少妇av免费看| 中国国产av一级| 国产精品99久久99久久久不卡 | 中文资源天堂在线| 日韩伦理黄色片| 精品少妇内射三级| 97在线视频观看| 伦理电影大哥的女人| 欧美激情国产日韩精品一区| 少妇人妻一区二区三区视频| 成人漫画全彩无遮挡| 国产片特级美女逼逼视频| 欧美3d第一页| 国产午夜精品一二区理论片| 免费大片黄手机在线观看| 我的老师免费观看完整版| 最新中文字幕久久久久| 精品少妇内射三级| 秋霞在线观看毛片| 精品人妻一区二区三区麻豆| 看免费成人av毛片| 制服丝袜香蕉在线| 26uuu在线亚洲综合色| 色哟哟·www| 国产精品99久久99久久久不卡 | 成人毛片60女人毛片免费| 午夜91福利影院| 国产精品熟女久久久久浪| 欧美日韩亚洲高清精品| 人人妻人人看人人澡| 欧美精品一区二区大全| 一本色道久久久久久精品综合| 亚洲成人手机| 特大巨黑吊av在线直播| 国产精品国产av在线观看| 午夜影院在线不卡| 美女中出高潮动态图| 亚洲三级黄色毛片| 日韩在线高清观看一区二区三区| 嘟嘟电影网在线观看| 嘟嘟电影网在线观看| 精品少妇黑人巨大在线播放| 麻豆成人av视频| 精品国产露脸久久av麻豆| 亚洲国产精品一区二区三区在线| 亚洲真实伦在线观看| 国产男女内射视频| 午夜老司机福利剧场| 国产美女午夜福利| 久久久久久久久久成人| 国内少妇人妻偷人精品xxx网站| 日韩大片免费观看网站| 2021少妇久久久久久久久久久| 国产美女午夜福利| 久久这里有精品视频免费| 黑丝袜美女国产一区| 久久人人爽人人爽人人片va| 日本av手机在线免费观看| 亚洲高清免费不卡视频| 热re99久久精品国产66热6| 日韩视频在线欧美| 成人亚洲欧美一区二区av| 成人毛片a级毛片在线播放| 久久精品国产自在天天线| 欧美少妇被猛烈插入视频| 老司机亚洲免费影院| 韩国av在线不卡| 能在线免费看毛片的网站| 少妇人妻精品综合一区二区| 乱人伦中国视频| 精品人妻熟女毛片av久久网站| 日韩精品有码人妻一区| 欧美精品国产亚洲| 精品久久久精品久久久| 久久久久久久国产电影| 日日摸夜夜添夜夜添av毛片| 亚洲av日韩在线播放| 日韩制服骚丝袜av| 女人精品久久久久毛片| 18禁在线无遮挡免费观看视频| 欧美3d第一页| 免费人成在线观看视频色| 日本黄大片高清| 国产亚洲欧美精品永久| 一区二区三区乱码不卡18| 国产黄色免费在线视频| 一边亲一边摸免费视频| 国产精品嫩草影院av在线观看| 黄色欧美视频在线观看| 成人毛片60女人毛片免费| 国产精品人妻久久久影院| 少妇精品久久久久久久| 欧美成人精品欧美一级黄| 亚洲精品乱码久久久v下载方式| 国产一区有黄有色的免费视频| 欧美日韩精品成人综合77777| 久久国产亚洲av麻豆专区| 国产欧美日韩一区二区三区在线 | 丰满人妻一区二区三区视频av| 色吧在线观看| 成人无遮挡网站| 免费黄色在线免费观看| 久久精品夜色国产| 日韩成人av中文字幕在线观看| 亚洲精品久久久久久婷婷小说| freevideosex欧美| 一级毛片 在线播放| 9色porny在线观看| 久久久久久久精品精品| 草草在线视频免费看| 国产精品成人在线| 人人妻人人澡人人看| 国产极品天堂在线| 欧美人与善性xxx| 久久婷婷青草| 高清视频免费观看一区二区| 国产欧美亚洲国产| 中国美白少妇内射xxxbb| 99热6这里只有精品| 欧美 亚洲 国产 日韩一| 国产免费一区二区三区四区乱码| 国产日韩一区二区三区精品不卡 | 午夜影院在线不卡| 久久精品国产自在天天线| 国产黄色免费在线视频| 女人久久www免费人成看片| 全区人妻精品视频| 免费av不卡在线播放| 久久热精品热| 男女边摸边吃奶| 一区在线观看完整版| 亚洲av成人精品一二三区| 国产永久视频网站| 国产一区二区在线观看日韩| 国产欧美另类精品又又久久亚洲欧美| 国产极品粉嫩免费观看在线 | 97在线视频观看| 人人妻人人澡人人看| 日韩一本色道免费dvd| 国产高清不卡午夜福利| 久久av网站| 亚洲精品一区蜜桃| 国产黄色视频一区二区在线观看| 亚洲四区av| 好男人视频免费观看在线| 久久午夜福利片| 亚洲精品第二区| av天堂久久9| 我的女老师完整版在线观看| 精品久久久久久久久av| 成人特级av手机在线观看| 在现免费观看毛片| 一本一本综合久久| 人人澡人人妻人| 一区二区三区四区激情视频| 国产日韩欧美视频二区| 久久久久久久久久久丰满| 三上悠亚av全集在线观看 | 精品人妻熟女av久视频| 国产在视频线精品| 亚洲欧洲日产国产| 国产成人精品婷婷| 国产熟女欧美一区二区| 97在线人人人人妻| 少妇精品久久久久久久| 成人午夜精彩视频在线观看| 国产深夜福利视频在线观看| 新久久久久国产一级毛片| 亚洲国产精品999| 岛国毛片在线播放| 国产av一区二区精品久久| 不卡视频在线观看欧美| 成人免费观看视频高清| 亚洲精品久久久久久婷婷小说| 最新中文字幕久久久久| 国产精品国产三级国产av玫瑰| 成人特级av手机在线观看| 精品久久久精品久久久| 亚洲国产最新在线播放| 免费观看无遮挡的男女| 亚洲四区av| 精品人妻熟女av久视频| 免费黄频网站在线观看国产| 青春草视频在线免费观看| 80岁老熟妇乱子伦牲交| 久久久国产精品麻豆| 日韩伦理黄色片| 噜噜噜噜噜久久久久久91| 精品少妇久久久久久888优播| 色婷婷久久久亚洲欧美| 国产欧美日韩一区二区三区在线 | 午夜老司机福利剧场| 国产中年淑女户外野战色| 黑人巨大精品欧美一区二区蜜桃 | 精华霜和精华液先用哪个| 亚洲精品久久午夜乱码| 国产一级毛片在线| kizo精华| 久久久久久久久久成人| 日日啪夜夜撸| 桃花免费在线播放| 嫩草影院新地址| 色视频www国产| 妹子高潮喷水视频| 国产一区有黄有色的免费视频| 国产爽快片一区二区三区| 国产高清三级在线| 久久精品国产亚洲av涩爱| 欧美精品人与动牲交sv欧美| 亚洲精品自拍成人| av在线观看视频网站免费| 国产一区亚洲一区在线观看| av播播在线观看一区| 午夜日本视频在线| 国产黄频视频在线观看| 亚洲美女视频黄频| 久久久久人妻精品一区果冻| 久久久国产一区二区| 亚洲怡红院男人天堂| 久久久久精品久久久久真实原创| 一级a做视频免费观看| 欧美丝袜亚洲另类| 免费少妇av软件| 中文天堂在线官网| 国产一区二区在线观看av| 国产黄片视频在线免费观看| 王馨瑶露胸无遮挡在线观看| 在线观看人妻少妇| 免费av不卡在线播放| 69精品国产乱码久久久| 美女xxoo啪啪120秒动态图| 亚洲精品国产av成人精品| 最近2019中文字幕mv第一页| 妹子高潮喷水视频| 国产精品麻豆人妻色哟哟久久| 日日摸夜夜添夜夜添av毛片| 国产白丝娇喘喷水9色精品| 国产精品久久久久久久电影| 亚洲精品亚洲一区二区| 久久久久视频综合| 久久国产亚洲av麻豆专区| 丝袜在线中文字幕| 精品国产一区二区三区久久久樱花| 久久99蜜桃精品久久| 国产成人免费观看mmmm| 能在线免费看毛片的网站| 日本黄色片子视频| 女人久久www免费人成看片| 亚洲国产精品专区欧美| 一级毛片我不卡| 精品少妇黑人巨大在线播放| 成人毛片60女人毛片免费| 一本—道久久a久久精品蜜桃钙片| 狂野欧美白嫩少妇大欣赏| 尾随美女入室| 欧美精品一区二区免费开放| 插阴视频在线观看视频| 高清毛片免费看| 日本黄大片高清| 免费久久久久久久精品成人欧美视频 | 黄色一级大片看看| 搡老乐熟女国产| 久久这里有精品视频免费| 看非洲黑人一级黄片| 丰满饥渴人妻一区二区三| 亚洲激情五月婷婷啪啪| 热re99久久国产66热| 国产伦精品一区二区三区四那| 国产精品三级大全| a级一级毛片免费在线观看| 欧美日韩视频高清一区二区三区二| 亚洲人成网站在线观看播放| 亚洲av成人精品一二三区| 日本黄色片子视频| 国产一区二区三区av在线| 男人狂女人下面高潮的视频| 国产极品天堂在线| 欧美xxxx性猛交bbbb| 街头女战士在线观看网站| 亚洲成人一二三区av| 国产免费福利视频在线观看| 亚洲av免费高清在线观看| 夜夜看夜夜爽夜夜摸| 中文字幕av电影在线播放| av福利片在线观看| 国产免费视频播放在线视频| 全区人妻精品视频| 日韩欧美精品免费久久| 日韩,欧美,国产一区二区三区| 一本久久精品| 国产一区二区三区综合在线观看 | 亚洲av日韩在线播放| 久久99热这里只频精品6学生| av视频免费观看在线观看| 久久久久久久精品精品| 国产在线免费精品| 色婷婷久久久亚洲欧美| 精品熟女少妇av免费看| 久久人人爽av亚洲精品天堂| 五月开心婷婷网| 精品少妇内射三级| 精品人妻偷拍中文字幕| 曰老女人黄片| 全区人妻精品视频| 五月天丁香电影| 亚洲欧美精品自产自拍| 国产中年淑女户外野战色| 欧美亚洲 丝袜 人妻 在线| 特大巨黑吊av在线直播| 成人黄色视频免费在线看| 亚洲精品成人av观看孕妇| 国产一区二区在线观看日韩| 精品人妻一区二区三区麻豆| 一本色道久久久久久精品综合| 久久精品国产亚洲网站| 中文字幕人妻丝袜制服| 曰老女人黄片| 大话2 男鬼变身卡| 亚洲精品第二区| 日日啪夜夜撸| 晚上一个人看的免费电影| 看非洲黑人一级黄片| 中文字幕人妻熟人妻熟丝袜美| 各种免费的搞黄视频| 丝瓜视频免费看黄片| 国产国拍精品亚洲av在线观看| 女的被弄到高潮叫床怎么办| 日韩av在线免费看完整版不卡| 亚洲av免费高清在线观看| 亚洲精品,欧美精品| 欧美xxxx性猛交bbbb| 日本午夜av视频| 久热久热在线精品观看| 欧美日韩一区二区视频在线观看视频在线| 日韩 亚洲 欧美在线| 免费高清在线观看视频在线观看| 黄色视频在线播放观看不卡| 亚洲天堂av无毛| 免费黄网站久久成人精品| 色吧在线观看| 国产色爽女视频免费观看| 啦啦啦在线观看免费高清www| av视频免费观看在线观看| 欧美性感艳星| 午夜福利在线观看免费完整高清在| av国产久精品久网站免费入址| 99国产精品免费福利视频| 最近中文字幕高清免费大全6| 亚洲久久久国产精品| 夜夜骑夜夜射夜夜干| 国产成人freesex在线| 综合色丁香网| 日本免费在线观看一区| 我要看黄色一级片免费的| 日韩人妻高清精品专区| 国产亚洲一区二区精品| 大话2 男鬼变身卡| 免费看日本二区| 69精品国产乱码久久久| 国产欧美日韩精品一区二区| 18禁在线播放成人免费| 国产精品一区二区在线观看99| 国产免费一区二区三区四区乱码| www.av在线官网国产| 国产免费一区二区三区四区乱码| 国产免费福利视频在线观看| 啦啦啦啦在线视频资源| 亚洲精品,欧美精品| 国产成人aa在线观看| 少妇被粗大猛烈的视频| 亚洲精品乱码久久久v下载方式| 成人影院久久| 亚洲欧美成人精品一区二区| 99精国产麻豆久久婷婷| 99久久精品国产国产毛片| 日产精品乱码卡一卡2卡三| 欧美精品亚洲一区二区| 欧美区成人在线视频| 男人爽女人下面视频在线观看| freevideosex欧美| 久久精品久久久久久久性| 国产一区二区在线观看av| 亚洲人成网站在线观看播放| 免费观看无遮挡的男女| 亚洲av二区三区四区| 最近的中文字幕免费完整| 亚洲欧美一区二区三区黑人 | 啦啦啦啦在线视频资源| 国产一区二区在线观看av| 一个人免费看片子| 亚洲一级一片aⅴ在线观看| 蜜桃在线观看..| 免费观看无遮挡的男女| 亚洲av国产av综合av卡| 国产精品国产三级国产专区5o| 男女啪啪激烈高潮av片| 最近的中文字幕免费完整| 2022亚洲国产成人精品| 91成人精品电影| 国产成人精品婷婷| av黄色大香蕉| 免费观看a级毛片全部| 亚洲精品乱码久久久v下载方式| 美女国产视频在线观看| 国产色爽女视频免费观看| 免费观看无遮挡的男女| 免费高清在线观看视频在线观看| 97在线人人人人妻| 啦啦啦中文免费视频观看日本| 国产探花极品一区二区| 在线 av 中文字幕| 亚洲av欧美aⅴ国产| av线在线观看网站| av又黄又爽大尺度在线免费看| 性高湖久久久久久久久免费观看| 日韩亚洲欧美综合|