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

    A comprehensive review of in-situ polymer gel simulation for conformance control

    2022-03-30 13:52:06BaojunBaiJianqiaoLengMingzhenWei
    Petroleum Science 2022年1期

    Baojun Bai,Jianqiao Leng,Mingzhen Wei

    Missouri University of Science and Technology,Rolla,65409,USA

    Keywords:In-situ gel Mechanism Simulation model Application studies

    ABSTRACT Gel treatment has been widely applied in mature oilfields to improve sweep efficiency and control water production.Correct numerical simulation is of major importance to the optimization design and prediction of a successful gel treatment.However,there exist many problems in current simulation studies in published literature.This paper first presents a comprehensive review on the major factors that have been considered at different gelation stages during gel treatment,the models used in the commercial/inhouse simulators,and current numerical simulation studies on both laboratory and field scales.Then we classify the current in-situ gel numerical simulation problems as 1,deficient model problem that has published numerical model but has not been applied in simulator and application studies;2,missing model problem that does not have published quantitative model;and 3,inaccurate application problem that does not consider the major factors of gel performance,based on the reasons from some questionable results of current simulation studies.Finally,we point out the major research efforts that should be made in the future to better simulate in-situ gel treatment process.The review indicates that numerous simulation studies using commercial software packages intend to predigest the gel treatment,many of which,however,ignore important mechanisms and mislead the operation of gel treatment.In fact,a full assessment of simulating in-situ gels cannot be achieved unless the quantitative models can be qualified in terms of transport and plugging mechanisms based on the experimental results.

    1.Introduction

    Excessive water production is a common problem as reservoirs mature due to reservoir heterogeneity and low sweep efficiency.For heterogeneous reservoirs,after thief zones(or water channels)are formed by extensive water floods,improved sweep is commonly needed immediately.Gel has been widely applied in oil industry including sealant for wellbore leakage (Zhu et al.,2021),fluid(Wu et al.,2021),and plugging agent for enhanced oil recovery(Zhao et al.,2021).Conformance control using gel treatment is a technique to encourage displacing movable oil in an un-swept zone and to improve water drive closer to optimal conformance condition (Bai et al.,2015).With proper application,gel treatment is proved to be an effective and economic solution to control the conformance and reduce water production from oil and gas fields.

    In-situ gel treatment was developed in 1970s (Needham et al.,1974) for conformance control.The concept of gel treatment is to improve the sweep efficiency of water flooding by placing a blocking slug of gel at water channel to seal the thief zone.Many types of gels have been tested in the past decades,the success experience favors the polymer-based gel system proposed by Sydansk (1988,1990) that used partially hydrolyzed polyacrylamides (HPAM) and Chromium III (Cr(III)) acetate for conformance control due to its relatively long and stable gelation time.Thus,this review will mainly focus on the HPAM/Cr(III)acetate type of gel.The basic operation for this technology is to pump the polymer solution(e.g.,0.5 wt%)and crosslinker system(e.g.,1/40 of polymer) as a mixture solution,called gelant,into the formation with a relatively low flow rate to reduce gelant invasion in matrix.Then,shut-in the well for certain time ensuring in-situ gelation taking place sufficiently.The gelation process refers to the transition from polymer/crosslinker to crosslinked polymer gel.During this period,the slug of three-dimension gel as a permeability modifier or barrier in the preferential water channel is formed.Finally,open the well and the water drive can be diverted to the unswept zone.Based on the location of gel placement,the treatment can be categorized by injection profile control,in-depth flow path diversion,and producer water shutoff.

    Based on experiences,Sydansk and Southwell (2000),Sydansk et al.(2005) and Seright et al.(2011a) stated that the conformance control using in-situ gel treatment is a very complicated concept that strongly depends on the gel-,reservoir-,well-,and formation-properties.On the other hand,based on their compositions and application conditions,effect of in-situ gel treatment also depends on the gel properties.To ensure a successful treatment and to optimize the performance of in-situ gel,numerous experiments have been carried out to study the mechanisms.However,due to limitation of budgets,instruments and research scale of lab experiments,the accurate critical conditions and optimization for gel treatments are commonly difficult to obtain.Thus,numerical simulation studies are vitally important to study the laboratory results and to optimize the field applications.However,due to complex physiochemical properties of polymer gel,quantitative studies of in-situ gel treatment are very challenging.Qualified simulation of gel treatment requires the consideration of major mechanisms,principal physio-chemical phenomena,and important influence factors;the consistency between laboratory and simulation results;and the capability of computational efficiency.

    Therefore,the purposes of this review are to summarize the mechanisms,to examine the capability of commercial/in-house simulators and to critically investigate the simulation studies.With this review,the researchers interested in gel treatment simulation can improve the simulators and simulation methods to meet current needs of mature oilfields.

    2.Mechanisms discussion of in-situ gel treatment for numerical simulation

    Simulation of in-situ gel is quite complicated because the system contains great change in viscosity and flow regime before and after gelation.Moreover,the system has both fluid and solid properties during crosslinking process and after gel is formed.It is necessary to categorize in-situ gel treatment as three stages based on the gelation process,including gelant stages,gelation stage and gel stage,separately.

    2.1.Gelant stage

    Gelant is a mixture of polymer,crosslinker and perhaps some additives.Before gelation,the gelant has quite the same flow behavior as polymer solution alone (Seright,1991a,1991b,1991b).The gelant properties and parameters to be considered for a numerical simulation include polymer rheology,polymer and crosslinker adsorption/retention,and inaccessible pore volume.

    Polymer rheology.Resistance factor(Fr)is a critical parameter to quantify the rheology of polymer in laboratory experiment.It is measured by pressure drop ratio of waterflooding to polymer injection at the same flow rate (Eq.(1)).Assuming no reduction of permeability caused by polymer,we can estimateFrby viscosity ratio due to the same effective permeability to gelant and water.If we assume water viscosity is unit,we can further estimateFrby the apparent viscosity of gelant.

    where ΔP is pressure drop,λ is mobility,kis effective permeability,μp,appis polymer apparent viscosity,subscript w refers to water and p refers to polymer.

    The rheology of polymer,as a non-Newtonian fluid,has been studied for decades by many researchers(Gogarty,1967;Graessley,1974;Gleasure,1990;Seright et al.,2008;Smith,1970;Zhang and Seright,2015).Seright et al.(2011b) summarized the rheology responses of polymers and concluded that the resistance factor could be a function of polymer concentration,salinity,shear rate(velocity),which depended on the type of polymer and formation properties.Their studies reported that Xanthan polymer only behaved a shear thinning rheology that the resistance factor decreases with increasing shear rate.HPAM polymer behaved a Newtonian flow regime at low shear rate,a shear thinning rheology at low to medium shear rate and a shear thickening rheology at medium to high shear rate.Shear thickening means that the resistance factor increases with increasing shear rate.Their studies also showed that the critical condition for the onsets of shear thinning and shear thickening depended on the polymer concentration,salinity,and porous media permeability.Besides,Seright et al.(2011b),Zechner et al.(2013)and Ma and McClure(2017)stated that if the formation contained open fractures,such as the natural fractures,extension of hydraulic fractures and wormholes,HPAM polymer would perform only shear thinning rheology,because the polymer elongations in small pores were not observed in large pores or fractures.

    Seright et al.(2011b) studied the impact of polymer concentration and salinity on polymer rheology.The results showed that with increasing polymer concentration from 200 to 1600 ppm,the onset of shear thickening(the velocity required to let polymer start to behave shear thickening) decreased by a factor of two.The results also showed that polymer tended to have a higher shear thickening coefficient (greater shear thickening slope) in high salinity water than in low salinity water.

    Polymer and crosslinker retention.During injection of gelant,polymer and crosslinker retention can have a major impact on gelant penetration into a reservoir and the gelant composition in different locations due to chromatographic effect.Adsorption was reported as the major mechanism for polymer and crosslinker retention that might occur via physical and chemical adsorption.The adsorption process can be characterized by determining the number of molecules adsorbed per surface area.μg polymer per g contaminant is the commonly used unit for polymer and crosslinker adsorption.Green and Whillhite (1998) reported the retention of polymer ranged from 9 to 700 μg/g for polymer concentration 500-3000 ppm,which greatly influenced the transport of polymer during placement.Langmuir adsorption isotherm is commonly suggested for the fitting of retention result of polymer.Polymer has a large molecule composed of repeating subunits bound together by covalent bonds.Thus,besides adsorption,Ferreira(2019),Zhang and Seright(2014)stated that retention data suggested the bridging,clogging and entrapment could also be the main mechanisms of polymer retention as a function of polymer hydraulic dynamic radius and pore throat size.The recent study (Dandekar et al.,2021) also indicated the polymer solution retention data was not fitted well with Langmuir adsorption isotherm.Stavland et al.(1993)studied the retention of crosslinker and stated that different retention mechanisms such as ion exchange and precipitation could be observed having alone or combination effect on crosslinker retention,which was related to pH value,carbonate concentration and core types.

    Inaccessible pore volume (IAPV).Dawson and Lantz (1972)observed the IAPV of polymer that some pore volume is accessible to small molecule water but inaccessible large molecule polymer.Since then several mechanisms conducting IAPV have been reported and summarized by Leng(2021):(1)Molecule sizes restricted polymer from penetrating into some pore space (Shah et al.,1978);(2) Existence of depletion effect of large molecules built up a depletion layer at pore walls and prevented the mass centers of polymer molecules from reaching pore walls(Chauveteau,1981;Omari et al.,1989;Sorbie,1989);(3)Unfavored entropic effect pushed polymers away from solid boundaries(Liauh et al.,1979);(4)In-situ rheology behavior of polymer was different from intrinsic rheology,which made polymer flow in-situ faster than expected (Chauveteau et al.,1984;Ferreira,2019;Stavland et al.,2010).Though the value and explanation of IAPV concept were not consensus,the existence of IAPV that influenced the transport of polymer was commonly acknowledged.Zhang and Seright (2014) summarized the IAPV data from the laboratory results reported in literatures and concluded that the IAPV value had a wide range from 0% to 48%.Their results showed that the IAPV was positively related to theSorand mole weight,negatively related to the resistance factor and retention.However,due to inconsistency of the IAPV results,no quantitative models were reported for IAPV model.

    2.2.Gelation stage

    The process of gelation contains inner crosslinking and outer crosslinking.Inner crosslinking refers to that single crosslinker connects itself to two adjacent polymer chains and forms two dimensional(2D)structures.Outer crosslinking refers to that these 2D crosslinker-polymer structures continue connecting themselves to adjacent 2D structures and form 3D structure gel.Commonly,the inner crosslinking takes much longer time than outer crosslinking.The crosslinking process is governed by gelation kinetics.Gelation kinetics has been studied by many researchers and many types of explanations have been proposed.Baylocq et al.(1998) explained the HPAM/Cr(III) acetate gelation by three steps’ reaction,which separate inner crosslinking into two steps including single crosslinker to one adjacent polymer chain and single to two adjacent polymer chains.They suggested that the triangular structured Cr(III)should hydrolyze in three steps.Jain et al.(2005)stated that the gelation process included two step-by-step reactions:uptake and crosslink.

    Gelation time is a fundamental parameter in oilfield applications that quantifies the total time for gel to form.However,due to different structures,the gelant resistance factor or dynamic viscosity will keep at low level at the stage of inner crosslinking but increase sharply at the stage of outer crosslinking.Thus,the transition from inner to outer crosslinking stage could be more important for field application and simulation studies.Sydansk(1988) used bottle test and defined gel-strength code to represent the stage of gelation.Winter and Mours (1997) applied Tung and Dynes method to define the starting time of outer crosslinking as the gel point.The method measured the elastic modulus (G′) and loss modulus (G′′) and used the time when the ratio ofG′′toG′is equal to 1 as the gel point.Romero-Zer′on et al.(2008)implemented nuclear magnetic resonance (NMR) to investigate the bulk relaxation rate during gelation and proposed the interception point of two stages of crosslinking could be more feasible and accurate to represent the gel point.

    Commonly,it takes a few hours for HPAM/Cr(III) type gel to reach gel point and another one or couple hours to reach gelation time(Sydansk et al.,2005).However,the recent publication stated the delayed technology of crosslinking so that the gel point could take up to one month (Cordova et al.,2008;Sun et al.,2016).

    The other key mechanism of gelation is the reaction rate of each crosslinking stage.Prud'homme et al.(1983) studied the HPAM/Cr(III) gel gelation process using rheological monitoring and proposed Arrhenius equation for the reaction rate.The reaction kinetics has polymer reaction order of 2.7 and crosslinker reaction order of 2 for the outer crosslinking stage.Scott et al.(1987)simplified the reaction kinetics by first order reaction for each reactant (second order in overall) so that the model could be generalized to more complex mixture of polymer and crosslinker.Romero-Zer′on et al.(2008) found that the reaction kinetics were not same for inner and outer crosslinking processes and the inner crosslinking was much slower than outer crosslinking.Their results proposed that only the outer crosslinking stage of gelation follows a second order overall reaction(first order on polymer and first order on crosslinker).

    2.3.Gel stage

    After a gel is formed,it is very difficult to transport through common porous media,but it might be able to propagate through fractures or fracture-like channels.The primary mechanisms in this process contain propagation and retention.Residual resistance factor is a key factor to evaluate the permeability reduction caused by retained gel.

    Propagation.Formed gel propagation through fractures or fracture-like channels has different transport behavior from gelant solution.Seright(1999,2001)studied bulk gel propagation through fractures and concluded that the pressure gradients required to propagate gels were greater than those for flow of gelants and a threshold pressure was required to mobilize the gel.

    Retention.Gel can retain during propagation through porous media or fractures.Gel retention may result from many mechanisms including size exclusion of large molecule,chemical adsorption or attachment by the rock surface,gravity segregation,diffusion to the pinch-out pores,and bridging (Chauveteau et al.,1998;You et al.,2013;Zhang and Seright,2014).Wang et al.(1981) considered a Langmuir adsorption model as a function of aqueous phase polymer concentration and salinity to quantify the retention concentration of gel.Seright (2009) studied the mechanism of gel retention and shown that additional to the monolayer adsorbed on rock surface,gel could aggregate at the pore throat and accumulate in the pore space,which indicated a pore filling behavior.Charoenwongsa et al.(2012) also suggested that the retention of in-situ gel should include both adsorption layer and solid entrapment layer.

    Dehydration.In fact,during gel propagation in fracture,the retention of formed gel is very complex.Formed gel can retain on surface of fracture to matrix.Moreover,due to superabsorbent property,the water molecules in gel structure can be squeezed out and flow through the fracture or to the matrix,which as a result increases the concentration of gel retained in fractures and forms a filter cake.This process is denoted as fluid leak-off in studies of bulk gel that has been reported by experimental researchers (Seright,1998,2001,2001;Brattek?s et al.,2018).

    Residual resistance factor.Due to gel retention in pore network,the permeability of the porous media can be reduced.Residual resistance factor (Frr) is commonly applied in laboratory to characterize the plugging efficiency,which is the permeability ratio before and after gel treatment shown in Eq.(2).

    wherekbis permeability before gel treatment andkais permeability after gel treatment.

    As measured in laboratory,Frris commonly a constant value that refers to the final result after gel is fully placed.However,in upscaled field application,the gel treatment requires the dynamic result with varied concentration and time.Thus,a function to relateFrrwith polymer concentration is necessary.Yuan (1991) studied the gel retention and used a tablet model to quantify the relationship between dynamicFrrand polymer concentration and the result was reported consistent with the exponential result.Stavland et al.(1994) proposed the permeability reduction as a linear function ofFrrand adsorbed polymer concentration.Cheng(2012) studied the retention of gel and stated that the retention of gel could cause higher permeability reduction than prediction of linear model.Their results employed an exponential function of retention concentration to fit the permeability reduction,which showed a better fitting.

    After gel is placed,several factors can reduce theFrrincluding chemical degradations,shear stress,and oil throughput.Seright(1988) studied the degradation of formed gel and found that unlike the erosion or desorption of polymer,after treatment of strong gel,the permeability reduction to water was stabilized rapidly and remained stable for more than six months.However,for long-term consideration,the chemical degradation that can gradually reduceFrrcannot be neglectable.Seright (2009) stated that shear degradation caused by the chasing water can be attributed by several aspects including the desorption of retained pore-filling gel due to increased pressure gradient.Brattek?s et al.(2016) applied visualization of bulk gel in fracture and observed the wormhole created in the gel slug due to higher shear stress in the center of fracture.Ganguly et al.(2002) applied coreflooding in tube model and reported that the rupture pressure gradient of gel structure was related to the tube size,length,salinity,gel composition and aging time.

    Besides,disproportionate permeability reduction (DPR) effect,or relative permeability modifier(RPM)effect,is also an important mechanism for gel treatment that most gel system can provide much higher permeability reduction to water than to oil.Thus,Frris different for water(Frrw)and for oil(Frro).Seright(2005)found that the oil throughput from 1 to 100 PV decreased the gel strength and increased the oil effective permeability gradually by factors from 5 to 10 times,which explained,as part of the reason,for the DPR effect.In fact,the proposed mechanisms for RPM are numerous.One of the most agreed is the segregated pathways for water and oil proposed by White et al.(1973) and Nilsson et al.(1998).

    3.Models for numerical simulation

    In this section,we will discuss the numerical models for in-situ gel simulation.The models include quantitative models that were published in literatures and models that were reported as build-in functions in commercial/in-house simulators including UTCHEM,SCORPIO,CMG STARS,POL-GEL,IORCoreSim and Eclipse 300.

    3.1.Gelant simulation

    Rheology:Polymer rheology is influenced by many factors including polymer concentration,shear rate,polymer plateau viscosity.The shear rate is influenced by flow velocity,rock phase porosity and permeability.Polymer plateau viscosity is influenced by polymer type and polymer concentration.

    In laboratory,polymer plateau viscosity is commonly measured using viscometer at 7.33 s-1.To fit the laboratory result,Thurston et al.(1987) proposed a linear regression function to fit the polymer plateau viscosity using polymer concentration and water viscosity (Eq.(3)).Stavland et al.(1994) added aFfactor that was a function of active polymer and crosslinked polymer.Delshad et al.(1996)added a salinity factor to consider the effect of effective salt concentration.

    CMG STARS definesF(x) function by Eq.(4) to modify the basic viscosity mixing rule and describes the change of viscosity that depends on mole fractionXof componenti.The modified model can quantify large scale non-linear viscosity change,which favors gel simulation.TheF(x)function is the input,and the simulator will calculate plateau viscosity using user-input.

    wherexis mole fraction;i,jare components in the water phase;pis polymer component;ncis the number of components in water.

    Polymer rheology can depend on the permeability and porosity of the porous media because the effective shear rate varies with pore size.To quantify the effective shear rate in porous media,Chauveteau and Zaitoun(1981)proposed shear rate as a function of u(1-φ)/(φk)0.5.Hirasaki and Pope (1974) proposed a function ofCannella et al.(1998) reported in two phase flow (oil and water)the effective shear rate(γe)as a function of shear coefficient,water phase velocity,water phase saturation,water phase effective permeability and porosity (Eq.(5)).Since then,though varied coefficients were reported,the general derivation of the Cannella model has not been changed and widely applied in each simulator.

    wherenis shear coefficient;Cis Cannella constant;uwis water velocity;Kis absolute permeability;krwis water relative permeability;Swwater saturation;φ is porosity.

    For rheology model,Camilleri et al.(1987) proposed Meter's equation (Meter and Bird,1964) to quantify the non-Newtonian behavior of polymer solution (Eq.(6)).Due to its simplicity and feasibility,the model is widely used in most simulators to simulate the shear thinning response of polymer rheology.

    where μappis apparent viscosity;γeqis equivalent shear rate;γsandPαare tunning parameters.

    Besides,Seright (1991a) summarized the rheology model of Xanthan and HPAM polymer including five shear thinning models which are power law (Bird et al.,1960),Carreau model (Bird et al.,1977),Chauveteau model(Chauveteau and Zaitoun,1981),Cannella model (Cannella et al.,1988),and Willhite empirical models(Willhite and Uhl,1988)and three shear thickening models which are Heemskerk dual power law model (Heemskerk et al.,1984),Hirasaki-Pope model (Hirasaki and Pope,1974),and Durst-Bird models (Durst et al.,1982).Delshad et al.(2008) combined the experimental results of Hirasaki and Pope(1974)and Masuda et al.(1992)with Carreau model and derived a dual power law model for HPAM polymer rheology (Eq.(7)).The model could quantify both shear thinning and shear thickening of polymer solution.The model also quantified the relationship between polymer viscosity at critical point (initial viscosity and μmaxmaximum viscosity)and polymer concentration,salinity.

    where μ∞is viscosity without polymer;μmaxis maximum viscosity of polymer solution under shear;γ is effective shear rate;λ,λ2,α are tunning parameters;τris shear stress;nis shear thinning coefficient andn2is shear thickening coefficient.

    For the IAPV model,the commercial simulators have their different aspects.Eclipse 300 (Eclipse 300 Manual,2014) assumes that the IAPV decreases the effective water saturation;POL-GEL(Yuan et al.,2000a) assumes that the IAPV reduces the polymer effective volume,which decreases the effective concentration of polymer;CMG STARS (Manual,2016) assumes that the IAPV decreases the effective porosity occupied by polymer.Although on different point of view,all simulators applied the constant value for IAPV of user input.

    For polymer and crosslinker retention,the most widely used is the Langmuir adsorption model that is a function of aqueous phase concentration and salinity (Langmuir,1918;Yuan et al.,2000a;Manual,2016) (Eq.(9)).

    Stavland et al.(1994) proposed more complex retention model for crosslinker that included ion exchange and precipitation (Eq.(10)).The model was later equipped in UTCHEM (Goudarzi,2015).

    3.2.Gelation simulation

    The simulation of gelation contains the transition from polymer/crosslinker to small aggregates(inner crosslinking)and to large 3D structural crosslinked gel (outer crosslinking).Scott et al.(1987)applied Arrhenius reaction rate (in Eq.(11)) to qua ntify the gelation kinetics.The method assumed 1 mol of polymer and 1 mol of crosslinker forms 1 mol of gel (m=n=1).With their method,the laboratory measured gelation time was converted to the reaction frequency (K) as an input and a function of temperature,activation function,gas constant and Arrhenius constant.The analytical solution provides the gelation time using Eq.(12).The method has been applied in CMG STARS and were used in many simulation studies.

    Stavland et al.(1994) considered this transition process as the changes in polymer properties that contained fraction of polymer following gel properties instead of using another component ‘gel’(Eq.(13)).Thus,in their model,the component number was reduced.

    whereFgis fraction of polymer that has changed to gel properties;Cis concentration;expp,expx,exphare tunning parameters;Eais activation function;Ris gas constant;T0is standard temperature andTis temperature in Kelvin.

    3.3.Gel simulation

    For gel propagation in fractures,Seright(1998)found that after a certain pressure was reached,theFrof formed gel kept constant which indicates a Bingham type fluid.The flow rate was calculated using Eq.(14).

    From one Shrove Tuesday to another, much may occur to weigh down the heart; it is the reckoning of a whole year; much may be forgotten, sins against heaven in word and thought, sins against our neighbor, and against our own conscience

    where q is flow rate;hfis fracture half length;wfis fracture width;μgis gel viscosity;Lis model length;Δp is pressure drop;y0is gel cake thickness.

    Wang and Seright (2006) applied a power law model to fit the relationship of flow rate andFrfor three different concentrations of polymer.The pressure drop was calculated using Eq.(15).

    wherea,n1,andn2are tunning parameters that are equal to 156,0.26 and 1.52 for 0.5%/0.0417% ratio of HPAM/Cr(III) gel.

    Ouyang et al.(2013) suggested Herschel-Bulkley rheology model for formed gel propagation with a better fitting result.The flow rate was calculated using Eq.(16).

    whereA=,τ0is yield stress;nisshear coefficient.

    For gel retention,Stavland et al.(1994)stated that gel retention could be quantified by Langmuir type equation (Eq.(17)).Their model considered a multiplier (Qm) that was related to fraction of polymer crosslinked.The model was claimed not fully correct in mechanisms due to particle clogging and entrapment but was reported capable of fitting several sets of experimental results.The model considered the constant rate of filtration.

    whereQm=Qp+CptFgsis aqueous fraction of gel;Ag0is minimum level ofFgs;AgaandBgaare tunning parameters for adsorbed gel and filtrated gel;Bg0represents the onset of more rigid gel formation(rapid increase in adsorption)when filtration occurs;is adsorption capacity;bis model input depend on salinity;Cp,aqis aqueous polymer;Cptis total amount of polymer present;φ,kare porosity and permeability,respectively.

    Unfortunately,no published simulation research studied the more sophisticate situation of retention.For reference,Lohne et al.(2010) studied the pore-throat trapping model and pore-lining retention model of drilling polymer fluid and derived Eq.(18) and Eq.(19),respectively,to quantify the trapping rateRtrapand lining retention rateRline(unit in PV).

    whereCtis polymer gel concentration at timet,σ is statistical result of trapping PV of polymer,uwis flow velocity,φ is porosity,Swis water saturation,λ1is trapping parameter,λ2is dragging parameter,λ3is adsorption parameter and λ4is desorption parameter,all parameters are in fraction of polymer retained per travel length unit.

    For fluid leak-off model,several simulation models are available.The conventional Carter's model (Eq.(20)) was proposed by Howard and Fast (1970).The model assumed the filter cake forms uniformly on the fracture surface when fluid leak-off.

    whereulis leak-off rate andtis leak-off time.

    Seright (2003) modified Carter's model and assumed that filter cake forms non-uniformly,which can better fit the experimental results.Brattekas et al.(2015) proved the Seright's model by the MRI observation of the gel cake formation.The leak-off rate is calculated using Eq.(21).

    where umis leak-off rate from mobile gel;wfis fracture width.

    Permeability reduction is contributed by the combination effect of the gel retention and dehydration.Based on the laboratory measuredFrr,the permeability reduction has been modeled as a mobility divider by many researchers.The most widely used model in publications is the linear function ofFrr(Eq.(22))(Manual,2016;Eclipse 300 Manual,2014;Lashgari,2018;Yuan et al.,2000a).The model only considered one phase permeability reduction.

    whereCadis adsorption concentration of gel;Cad,mis adsorption capacity.

    Stavland et al.(1994) has proposed a more complex permeability reduction model that considered porosity,permeability and the fraction of gel formed (Fg) (Eq.(23)).

    whereaFrr,bFrrare tunning parameters;Cg,ads,Cp,adsare adsorption concentration of gel and polymer.

    IORCoreSim simulator considers the permeability reduction as a function of effective water porosity and IAPV and the shear factor.The equation is shown in Eq.(24).

    wherefshis shear factor that is a function of shear rate based on Carreau-Yasuda model;frkwis tunning factor;φswis swelling factor of gel.

    For degradation of formed gel as a plugging agent,the mechanism contains chemical degradation and shear degradation.The chemical degradation rate (Rcd) is commonly simulated using first order(m=1)Arrhenius type kinetics(Eq.(25)).Due to low rate of chemical degradation,this model is capable of simulating the time dependentFrrreduction.

    whereKcdis reaction frequency of chemical degradation.

    The shear degradation has been studied by several researchers.Scott et al.(1987) considered the shear degradation as the mass reduction of formed gel using a constant rate for user input.Stavland et al.(1994)considered the effect of shear degradation by a multiplier onFrr(Frrf) (Eq.(26)).

    where Φ0is minimum pressure gradient to initiate shear degradation;Φ is effective pressure gradient;and Φfis tunning parameter.

    Lohne et al.(2017) considered the degradation as the mole weight reduction that was a function of effective shear stress (τ),degradation rate (rdeg),current mole weight (Mg),polymer hydraulic radius (Rp) (Eq.(27)).

    where αdis the tuning parameter.

    As a summary,Table 1 show the models that should be considered during gel treatment process.The major mechanisms refer to macroscopic mechanisms that commonly can be observed in laboratory as phenomena.The minor mechanisms refer to mechanisms behind these phenomena that need interpretation in quantitative models.All the simulation models are listed in column four and five,with considered factors and corresponding mechanisms.The ‘NA’ means no models are reported in the simulators.The last column describes the eligibility of the models in commercial simulators.The eligible models marked by ‘Y’ mean the simulation models in current simulators can provide effective estimation of the corresponding mechanisms.On the other hand,the ineligible models marked by‘N’mean either the models are not available in simulators,or the quantitative models have not been studied.

    Table 1 Summary of mechanisms,simulation models and concerned factors for in-situ gel simulation.

    Table 1 (continued)

    Table 2 Application studies review of in-situ gel simulation.

    Table 2 (continued)

    4.Application studies of in-situ gel simulation

    We summarized the published application studies of in-situ gel simulation as shown in Table 2.The first three literatures conducted the sensitivity studies on gelation kinetics and gelant injection dynamics using simulation models in laboratory scale.The influencing factors included pH value,temperature,polymer/crosslinker concentration,residence time,gelation time,Frrand gel volume.The gelation process determines the amount of gel that has been formed and retained.The gelant viscosity and rock permeability can be influenced by the fraction of formed gel.The process is quite complicated and is very difficult to be examined using laboratory experiments but has been studied using simulation.Helleren(2011) studied the gelant injectionFr,gel retention and post-flushFrras a function of the factors influencing gelation kinetics.If the gel formed during the injection of gelant,theFrrdue to gel retention will greatly influence the dynamics of gelant injection.Hatzignatiou et al.(2014)studied this effect using a simulation core model.They studied the injectionFrusing constant flow rate and injection rate variation using constant pressure,respectively.The gelation time also influences the loss of gel and further the placement depth.Hadi Mosleh et al.(2016) studied the effect of gelant composition and pH value on gelation kinetics and its effect on the gelant penetration.

    The other literatures listed in Table 2 stated field scale applications.For objective function,cumulative oil recovery is considered in every simulation study as the key result of gel treatment.However,as a short-term conformance control method,in-situ gel treatment may not influence cumulative oil recovery as much as other method such as polymer flooding.Thus,the result of cumulative oil recovery may reduce the significance of in-situ gel treatment.Compared with cumulative oil recovery,oil rate and water cut may reflect the short-term result of gel treatment more accurately.Most of the literatures included water cut as one of the major results but only Herbas et al.(2004a,b)considered the oil rate as an objective function to evaluate gel treatment.On the other hand,water(injection/production)profile is a critical and concise method to evaluate the efficiency of in-situ gel treatment.Unfortunately,only three literatures including Hughes et al.(1990),Hwan (1993)and Menzies et al.(1999) considered the profile change due to gel treatment.For injection operation,Hughes et al.(1990) and Seright et al.(2012)considered the zonal isolation injection,while the others applied ‘bullhead’ injection of gelants.

    For influencing factors,many factors have been discussed in the sensitivity analysis parts of these application studies.However,the methods to apply these factors might not reflect the mechanisms that have been discussed earlier in this paper.For example,Frrwas considered in all the literatures,but theFrrwas only related to retained gel concentration.Moreover,theFrrdistribution results showed that most of the literatures reported very lowFrrvalue in each grid (e.g.,lower than 100) after treatment,though some claimed maximumFrrover 104as input.As one of the explanations,Lee and Lee(2013)discussed the significance of using refined local grids and heterogeneities at treatment zone especially for the nearwellbore treatment and concluded that theFrrcould be unexpected low due to dilution of formation water without local grids refinement.Another example is the consideration of DPR effect.Both Sheshdeh et al.(2016)and Alfarge et al.(2018)considered the effect of DPR in simulation,however,the methods they applied did not consider the practicalFrromeasured from labs.Alfarge et al.(2018)investigated the RPM effect on gel plugging efficiency by using different permeability reduction to water (RKW)/permeability reduction to oil (RKO) ratio.However,the RKO was assumed to be one for the study which underestimated the formation damage caused by retained gel.Sheshdeh et al.(2016)applied permeability reduction by an interpolation method that two different relative permeability curves were set for water flooding before gel placement and after gel placement,respectively.This method could quantify the different permeability reduction to water and oil by retained gel after placement by assuming the constant permeability reduction in each grid.Thus,it inevitably ignored the transient period between two sets of relative permeability curves.This is critical because the permeability reduction in each grid was also a function of gel concentration.

    5.Discussion and suggestions

    As reviewed above,the simulation studies of in-situ gel treatment have been conducted for decades.Plenty of models have been proposed by researchers to quantify the transport and plugging performance of in-situ gel from gelant to formed gel.Based on these fundamental studies,numerous sensitivity analysis and optimization articles have been reported for the guidance of field scale application of in-situ gel treatment (Gao et al.,1993;Ghahfarokhi,2016;Goudarzi,2015;Herbas et al.,2004a,b;Khamees and Flori,2018;Khamees,2020;Lee and Lee,2013;Temizel et al.,2016;Yuan et al.,2000b).However,all of the simulations provided optimistic results,more or less,in terms of oil recovery and water cut.Meanwhile,several theoretical studies suggested the fastidious reservoir environment required for effective gel treatment (Sydansk and Seright,2007).Seright et al.(2011a) stated the challenging operations for in-situ gel treatment and limited situations that in-situ gel treatment could be favored than polymer flooding.Besides,the field applications also not always reported the successful results (Aldhaheri et al.,2016).

    Based on the comprehensive review,we found the results of many numerical simulation studies were not quite convincible.As shown in Table 1,many key mechanisms do not have proper models in commercial simulators,such as IAPV of polymer flow,gel point in gelation stage,gel transport in porous media,gel retention,gel propagation,dehydration,and DPR effect in gel stage,etc.Some of these mechanisms have numerical simulation models published but have not been applied in simulators,the others do not have eligible quantitative models for numerical simulation.Besides,some mechanisms are not considered in application studies.Thus,we categorize the reasons for questionable results of application studies as the deficient model problems,the missing model problems,and the application problems.

    Deficient model:Deficient models refer to those that related mechanisms have been considered and are available in current commercial/in-house simulators;however,the models are not eligible to quantify the mechanisms properly or thoroughly.These models include gelation process,gel retention,andFrrmodels.

    The gelation process is considered in simulators using Arrhenius reaction model.However,the model of gel point has never been considered in the gelation process,which may greatly influence the gel placement.As mentioned in gelation stages section,the gel point marks the time the outer reaction takes place and the resistance factor starts to rise significantly.For current gelation models in simulator (Polymer+Crosslinker >Gel),the main problem is that once gelant is injected,the gel starts to form at highest reaction rate based on Arrhenius reaction kinetics.Simultaneously the early injected gelant can form gel and retain in the pathway.Therefore,the later injected gel will be diverted from the gel pathway (most likely the high permeability water channel).Consequently,the early injected and retained gel inevitably would increase the injection pressure due to permeability reduction.As a result,without consideration of gel point,the gel placement cannot be properly estimated,and the injection pressure cannot be properly matched.Many researchers discussed methods to determine the gel point;however,no publication indicates the numerical modeling of gel point.

    For retention of formed gel,the conventional Langmuir adsorption as shown in Eq.(17) was reported not accurate by several researchers (Seright,2000;Cheng,2012;Charoenwongsa et al.,2012).This is because beside the adsorption,gel can be trapped by pore throat,retained,and fill the pore.Therefore,the retention volume and plugging efficiency can be higher than estimated using monolayer Langmuir model and linear retention model in commercial simulators.Many simulation models were proposed;however,the models were derived theoretically and only considered single or limited mechanisms.Thus,no retention models could history match the amount of gel retained considering the adsorption,solid entrapment,and dehydration mechanisms.In field of polymeric drilling fluid,the retention model proposed by Lohne et al.(2010) may be partially qualified in terms of factors concerned.However,the application of their model still needs consideration of dehydration and needs more experimental results for validation.

    Correspondingly,theFrrmodel as a function of retention is also not eligible.Without consideration of gel entrapment and dehydration effect on gel retention,the simulation model ofFrrmay very possibly overestimate the damage in matrix and underestimate the plugging efficiency in channels.The general model to quantifyFrrin simulation studies is Eq.(22),which shows a linear relationship with retention.In fact,theFrris very sensitive to the gel strength and formation properties,such as permeability and porosity.Weak gel restricts flow in low permeability rocks by a factor that is the same or greater than that in high permeability rock;while strong gel reduces permeability of all rocks to the same low value (e.g.,micro-Darcy level) (Seright,2002).Thus,Frrmodel should be a function of gel strength,retention,velocity,channel types and permeability,and the other formation properties.

    Missing model.Some mechanisms do not have quantitative models such as IAPV model of gelant flow and formed gel transport models.As stated in mechanism section,IAPV is influenced by many factors.The constant input setting in the commercial simulators cannot properly estimate the effect of IAPV on gel penetration in low permeability strata.Thus,the simulation result may underestimate the damage of gel in oil bearing zone.However,though many researchers reported the IAPV effect,the studies were only qualitative.Based on our review,only one article discussed the derivation of IAPV that was reported by Ferreira and Moreno(2019).Their equation was based on the Stavland et al.(2010)'s rheology model with some assumptions.However,the model can hardly simulate the effect of IAPV on polymer flow because the model only fit the situation of shear thinning and assumes low viscosity of polymer that is close to water (RF~1).

    When gel is formed,researchers only discussed the transport regime in open fractures.The gel is commonly considered not capable of penetrating into the porous media with complex pore networks (Bai et al.,2015),so that gel will selectively flow to the fractures or fracture-like channels.However,in relative high permeability matrix,such as high permeable unconsolidated-sand porous media,weak gel might still be movable when the threshold pressure gradient is achieved.Thus,more research is necessary for transport regimes and models of formed gel in relative high permeability porous media and critical conditions for selective penetration as a function of pore size,velocity/shear rate,and gel strength.For gel propagation in fractures,several numerical simulation models have been proposed in literatures,but the simulators do not have capable models to quantify the gel transport.In high permeability fractures or fracture-like channels,formed gel can propagate with a much higher resistance factor and dehydrate on fracture-matrix contact to form a filter cake.The gel cake can decrease the crossflow of mobile gel from fracture to matrix.Therefore,without consideration of these transport mechanisms,the simulation result may erroneously estimate the placement of gel in terms of gel consumption and penetration.

    As discussed above and listed in Table 1,the problems of simulators include the ineligibility of models for gelation kinetics,gel retention andFrr;and the lack of models for IAPV and gel transport.

    Problems and suggestions for application studies of in-situ gel simulation.Based on the previous review and discusses,it seems that the simulation models in simulators most likely have underestimated the effect of gel for conformance control.The question is why the application studies of in-situ gel simulation always provided more optimistic results than analytical and field results,but the results were not very convincible.We concluded several problems that widely exist in the application studies.

    One of the problems is the lack of correct consideration of DPR effect.As discussed in mechanism section,many polymer gels can reduce permeability to water more than that to oil or gas (Alfarge et al.,2018;Liang et al.,1995;Seright,2009;Willhite et al.,2002);however,the permeability reduction to oil or gas cannot be neglected.Seright (2009) summarized theFrroresults from experiments and found that theFrroranged from 2.7 to 59 in Berea sandstone.Though theFrrois smaller thanFrrw,considering the irreversible retention of formed gel in porous media,the damage to oil bearing zone,withFrroup to 59,is severe.

    Theoretically,DPR effect was reported only of value for gel treatment at production well for water shutoff (Sydansk and Seright,2007),because only for production well treatment,the remaining oil would flow through the gel zone.However,in simulation studies,the DPR effect might also have influenced the mobile oil flow for in-depth treatment.This is because numerous literature applied excessive large volume of gel treatment(over one month gelant injection)for in-depth fluid diversion,which caused the gelant flowing in-depth through channels and invading the oilbearing matrices from channels.In these simulation studies,Frrowere ignored.Thus,even though a large amount of gel had been erroneously placed in oil-bearing strata,the DPR effect on oil effective permeability was not considered and damage was neglected.Thus,the simulation results were not convincible.

    Another problem is the lack of consideration of the matrix fracking or fracture propagation(or extension).Due to pressure rise caused by low mobility of gelant and gel retention,the injection pressure required for gelant placement can easily overcome the formation fracking pressure and cause the creation of new fractures in matrices or extension of existing fractures to the matrices(Khodaverdian et al.,2009).Consequently,the simulation of gel treatment needs consideration of geomechanics.With our review,we found that in most current simulation studies,the injection pressure was commonly set to be unlimited to ensure a good match of injection rate.Without consideration of fracking model during gelant placement,the injection pressure will increase much higher than practical result in field and cannot guide applications.Moreover,the simulated injection pressure will very likely be higher than pump limitation.Thus,although the simulation studies showed optimistic oil recovery and water cut reduction,the results could hardly provide reasonable estimations to the field operations.

    The third problem is the misconception of gelation kinetics.As discussed in previous section,the gelation contains an inner crosslinking (low viscosity period) and an outer crosslinking (high viscosity period).Current application studies did not consider the gel point and applied an Arrhenius equation to quantify the gelation process.The problem with this model is that highest reaction rate occurs at the beginning because the rate is positively dependent on the remaining reactants’ concentration.Thus,the gel will be quickly formed from the beginning of gelant injection,and the fluid viscosity will be increased.When a large amount of gelant was injected,in most application studies,the gelant were erroneously like a super-polymer-flooding agent instead of a plugging agent,which would greatly mislead the field operations (Seright and Brattekas,2021).Because improper low values ofFrrwere applied in these application studies,the oil recovery and water cut results seemed very optimistic.However,if retention model andFrrwere considered properly,the later injected gelant would be diverted to the low permeability matrices and damaged the formation due to early retention of formed gel.Because in laboratory experiments,the injection time is commonly shorter than the gel point,the measuredFrresult does not includeFrr,however,in field application,the injection time could be much longer.Thus,for gelation kinetics,we suggest applying simulations case-by-case considering field scale gelant placement time.For placement duration much shorter than gel point(small treatment amount),Frrmodel should not be included in gelant placement stage.Frmodel could be applied for rheology model of gelant directly andFrrmodel should be applied for permeability reduction model for post flsh,respectively.For gelant placement duration close to or smaller than gel point (obvious viscosity rise and retention will happen during gelant injection),bothFrmodel andFrrmodel should be applied in gelant placement andFrrmodel should be applied in post flush as well.

    Last but not least,for severe channeling (e.g.,superpermeability channel,fracture-like channel,or fractures),the fluid flow in channels should not be considered as Darcy flow.The high flow velocity in severe channels commonly exceed the Reynolds number(Re)for laminar or Darcy flow.Non-Darcy behavior is important for describing fluid flow in situations where super-high velocity occurs.Although the criticalRefor non-Darcy flow is very inconsistent that ranges from 0.1 to 1000 depending on the rock and fluid properties(Chilton and Colburn,1931;Ergun and Orning,1949;Fancher and Lewis,1933;Ma and Ruth,1993),a typical value ofRe=10 may be applied for gelant flow in severe channels.Hassanizadeh and Gray (1987) explained this value for non-Darcy flow as the result from the increase in the microscopic viscous force at high velocity.We suggest future studies to consider the non-Darcy flow in simulation of gelant flow in severe channels and implement appropriate models in simulators.

    6.Conclusion

    Simulation of in-situ gel treatment is a comprehensive work due to the complex physiochemical properties of gel and the interactions between gel and formation.In this review,we classified the process of in-situ gel treatment as three stages and summarized the major factors to be considered during each stage,the mechanisms for in-situ gel treatment and the published numerical simulation models.

    A comprehensive table of simulation models is provided to compare the essential mechanisms with the published simulation models and the eligibility in simulators.The results show that the mechanisms including plateau viscosity,effective shear rate,rheology,polymer/crosslinker retention,gelation time,and chemical degradation of formed gel,have eligible models in integrated commercial/in-house simulators.The mechanisms including gelation kinetics,gel retention andFrrmodel have corresponding models in simulators,but the models cannot fulfill the needs required by the mechanisms.The mechanisms including IAPV of gelant,formed gel transport model,do not have eligible models in simulators.The gel propagation in fractures has several published empirical models (as shown in Table 1) but has not been incorporated in simulators,while the critical condition for gel selective penetration still needs more studies on quantitative models.

    A critical review of the published application studies using commercial/in-house simulators shows that the results of these application studies are questionable because(a)the simulators lack the essential models for in-situ gel treatment and(b)the problems existing in simulations include incorrect consideration of DPR effect for in-depth fluid flow diversion treatment and water shutoff treatment,the lack of consideration on the matrix fracking or fracture propagation,the misconception of gelation kinetics,and the non-Darcy flow in severe channels.

    Simulation of in-situ gel has been studied by many researchers for decades,however,major improvements are still in need to properly quantify the process of gelant placement,gelation kinetics,and formed gel mechanisms for conformance control.

    午夜久久久久精精品| 每晚都被弄得嗷嗷叫到高潮| 亚洲乱码一区二区免费版| 免费在线观看完整版高清| 性色av乱码一区二区三区2| 日韩有码中文字幕| 99国产精品99久久久久| 久久精品91蜜桃| 日本 av在线| 亚洲性夜色夜夜综合| 操出白浆在线播放| 国产激情偷乱视频一区二区| 此物有八面人人有两片| 亚洲中文av在线| 给我免费播放毛片高清在线观看| 国产精品99久久99久久久不卡| 每晚都被弄得嗷嗷叫到高潮| 99热只有精品国产| 又大又爽又粗| 最新美女视频免费是黄的| 两性午夜刺激爽爽歪歪视频在线观看 | 少妇人妻一区二区三区视频| 久久香蕉国产精品| 1024视频免费在线观看| 亚洲天堂国产精品一区在线| 在线观看免费视频日本深夜| 亚洲男人天堂网一区| 午夜福利在线观看吧| 狂野欧美白嫩少妇大欣赏| 黄色丝袜av网址大全| 亚洲一区中文字幕在线| 婷婷六月久久综合丁香| 男人的好看免费观看在线视频 | 91麻豆精品激情在线观看国产| 亚洲av熟女| av片东京热男人的天堂| 韩国av一区二区三区四区| 97碰自拍视频| 天天添夜夜摸| 精品欧美一区二区三区在线| 成年版毛片免费区| 亚洲成人国产一区在线观看| 成人手机av| xxx96com| 午夜福利高清视频| 国产一区二区激情短视频| 1024手机看黄色片| 欧美高清成人免费视频www| 久久久久久人人人人人| 欧美性猛交╳xxx乱大交人| 国产成人一区二区三区免费视频网站| 免费在线观看日本一区| 欧美三级亚洲精品| 五月玫瑰六月丁香| 国产高清视频在线观看网站| 国产区一区二久久| 18禁裸乳无遮挡免费网站照片| 18禁裸乳无遮挡免费网站照片| 美女 人体艺术 gogo| 成年女人毛片免费观看观看9| 日韩欧美国产在线观看| 欧美日本视频| 国产熟女xx| 国产精品乱码一区二三区的特点| 在线观看www视频免费| 亚洲自拍偷在线| 亚洲,欧美精品.| 亚洲精品粉嫩美女一区| 成人午夜高清在线视频| 变态另类成人亚洲欧美熟女| 无限看片的www在线观看| 久久香蕉精品热| 亚洲精品一区av在线观看| 草草在线视频免费看| 黄色视频,在线免费观看| 男人的好看免费观看在线视频 | 在线国产一区二区在线| 人妻丰满熟妇av一区二区三区| 国产精品av视频在线免费观看| 亚洲av成人不卡在线观看播放网| 免费av毛片视频| 在线免费观看的www视频| 国产欧美日韩精品亚洲av| 91av网站免费观看| 久久久久久久精品吃奶| 1024手机看黄色片| 校园春色视频在线观看| 色综合婷婷激情| 亚洲精品在线美女| 一级a爱片免费观看的视频| 亚洲,欧美精品.| 亚洲精品中文字幕一二三四区| 精品欧美国产一区二区三| 亚洲精品国产一区二区精华液| 亚洲中文日韩欧美视频| 国产精品1区2区在线观看.| 免费在线观看影片大全网站| 精华霜和精华液先用哪个| 十八禁网站免费在线| 夜夜爽天天搞| 97超级碰碰碰精品色视频在线观看| xxx96com| 久久国产精品影院| 日韩成人在线观看一区二区三区| 国产精品久久久久久人妻精品电影| 亚洲国产中文字幕在线视频| 亚洲熟女毛片儿| 88av欧美| 一进一出抽搐动态| 岛国在线观看网站| www.999成人在线观看| 女警被强在线播放| 99久久99久久久精品蜜桃| 成人三级黄色视频| 青草久久国产| 欧美日韩国产亚洲二区| 他把我摸到了高潮在线观看| 国产成人av教育| 久久久久免费精品人妻一区二区| 国产精品影院久久| 两个人视频免费观看高清| 精品久久久久久成人av| 亚洲国产高清在线一区二区三| 国产av不卡久久| 国产日本99.免费观看| 搞女人的毛片| 丰满人妻一区二区三区视频av | 亚洲全国av大片| 久久精品综合一区二区三区| 欧美中文日本在线观看视频| 久久这里只有精品19| 50天的宝宝边吃奶边哭怎么回事| 97碰自拍视频| 久久亚洲真实| 日日爽夜夜爽网站| 少妇人妻一区二区三区视频| 俄罗斯特黄特色一大片| 精品久久久久久久末码| 啪啪无遮挡十八禁网站| 国产一区二区三区视频了| 18禁观看日本| 久久久久亚洲av毛片大全| 亚洲,欧美精品.| 一级片免费观看大全| 欧美色欧美亚洲另类二区| 99精品在免费线老司机午夜| 精品不卡国产一区二区三区| 丁香欧美五月| 亚洲九九香蕉| 99久久国产精品久久久| 脱女人内裤的视频| 亚洲精品中文字幕在线视频| 欧美黑人欧美精品刺激| 国产精品久久视频播放| 九色国产91popny在线| 身体一侧抽搐| 亚洲av熟女| xxx96com| 啦啦啦免费观看视频1| 欧美成人一区二区免费高清观看 | 午夜日韩欧美国产| 无遮挡黄片免费观看| 性色av乱码一区二区三区2| 国产久久久一区二区三区| 一本综合久久免费| 婷婷六月久久综合丁香| 麻豆成人午夜福利视频| 婷婷亚洲欧美| 精品久久久久久久毛片微露脸| 成人亚洲精品av一区二区| 精品国产超薄肉色丝袜足j| 天天添夜夜摸| 精品一区二区三区av网在线观看| 亚洲精品在线观看二区| 丝袜人妻中文字幕| 欧美 亚洲 国产 日韩一| 婷婷六月久久综合丁香| 97超级碰碰碰精品色视频在线观看| 欧美精品亚洲一区二区| 欧美性猛交黑人性爽| 精品久久蜜臀av无| 欧美三级亚洲精品| 日日干狠狠操夜夜爽| 久久国产精品人妻蜜桃| 国产精品1区2区在线观看.| 亚洲美女视频黄频| 欧美黑人精品巨大| 国产精品亚洲美女久久久| 亚洲va日本ⅴa欧美va伊人久久| 午夜福利高清视频| 19禁男女啪啪无遮挡网站| tocl精华| 淫妇啪啪啪对白视频| 我的老师免费观看完整版| 亚洲成人国产一区在线观看| 色精品久久人妻99蜜桃| 精品欧美国产一区二区三| 两个人看的免费小视频| 欧美三级亚洲精品| 亚洲精品色激情综合| 18禁美女被吸乳视频| 黄色片一级片一级黄色片| 午夜精品一区二区三区免费看| 无限看片的www在线观看| 亚洲精品中文字幕在线视频| 国产精品久久久久久精品电影| a级毛片在线看网站| 久久久久国产一级毛片高清牌| 黄色毛片三级朝国网站| 精品少妇一区二区三区视频日本电影| 日韩精品免费视频一区二区三区| 99在线人妻在线中文字幕| 日韩高清综合在线| 久久午夜综合久久蜜桃| 在线观看免费视频日本深夜| 国产一区在线观看成人免费| 欧美在线一区亚洲| 国产激情久久老熟女| 在线观看66精品国产| 中文字幕人成人乱码亚洲影| 99国产精品一区二区蜜桃av| 老司机福利观看| 黄色成人免费大全| 欧美性猛交╳xxx乱大交人| 精品久久久久久成人av| 成年版毛片免费区| 一区二区三区高清视频在线| www.999成人在线观看| 这个男人来自地球电影免费观看| 国产成人啪精品午夜网站| 全区人妻精品视频| www.999成人在线观看| 又大又爽又粗| 国产av麻豆久久久久久久| 国产午夜福利久久久久久| 老汉色av国产亚洲站长工具| av视频在线观看入口| 日韩欧美在线二视频| 国产精品1区2区在线观看.| 午夜福利欧美成人| 亚洲精品久久国产高清桃花| 成人三级做爰电影| 欧美不卡视频在线免费观看 | 欧美乱妇无乱码| 淫妇啪啪啪对白视频| 国产黄片美女视频| 19禁男女啪啪无遮挡网站| 精品一区二区三区四区五区乱码| 国产精品一区二区三区四区免费观看 | 两个人免费观看高清视频| 麻豆av在线久日| 欧美另类亚洲清纯唯美| 国产精品1区2区在线观看.| 国产成人系列免费观看| 俄罗斯特黄特色一大片| 国产高清激情床上av| 亚洲人成77777在线视频| 国产视频内射| 国内精品一区二区在线观看| 美女午夜性视频免费| 国产精品久久视频播放| 日本一本二区三区精品| 精品电影一区二区在线| 黄色视频,在线免费观看| 熟女少妇亚洲综合色aaa.| 国产亚洲精品久久久久久毛片| 岛国视频午夜一区免费看| 香蕉国产在线看| 国内精品久久久久精免费| 国产一区二区三区视频了| 久久久久久大精品| 中出人妻视频一区二区| 天天一区二区日本电影三级| 国产精品香港三级国产av潘金莲| 亚洲五月天丁香| 99热这里只有精品一区 | 欧美久久黑人一区二区| 欧美日韩精品网址| 久久精品夜夜夜夜夜久久蜜豆 | 欧美3d第一页| 国产亚洲精品一区二区www| 午夜福利欧美成人| 高清毛片免费观看视频网站| 黄色a级毛片大全视频| 亚洲 国产 在线| 免费看美女性在线毛片视频| 亚洲avbb在线观看| 欧美三级亚洲精品| 欧美性猛交黑人性爽| 老司机福利观看| 亚洲国产高清在线一区二区三| 九色成人免费人妻av| 欧美黄色淫秽网站| 变态另类成人亚洲欧美熟女| 99热只有精品国产| 中文字幕精品亚洲无线码一区| 国产免费男女视频| 丰满的人妻完整版| 可以在线观看的亚洲视频| 在线十欧美十亚洲十日本专区| 亚洲熟妇熟女久久| 亚洲精品国产一区二区精华液| 久久久久国产精品人妻aⅴ院| 亚洲无线在线观看| 日日爽夜夜爽网站| 国产一区二区激情短视频| 久久久久久久久免费视频了| 日韩欧美三级三区| 欧美久久黑人一区二区| 欧美黄色片欧美黄色片| 国产av一区在线观看免费| 亚洲国产精品sss在线观看| 亚洲人成伊人成综合网2020| 少妇粗大呻吟视频| 国产亚洲av嫩草精品影院| 久热爱精品视频在线9| 搡老妇女老女人老熟妇| 99国产精品一区二区三区| 亚洲性夜色夜夜综合| 精品久久久久久久人妻蜜臀av| 免费在线观看完整版高清| 黄色丝袜av网址大全| av欧美777| 亚洲欧美精品综合久久99| 无遮挡黄片免费观看| 久久久久性生活片| av中文乱码字幕在线| 成人av在线播放网站| 亚洲 欧美 日韩 在线 免费| 日韩高清综合在线| 夜夜看夜夜爽夜夜摸| 欧美黑人巨大hd| 两性午夜刺激爽爽歪歪视频在线观看 | 99riav亚洲国产免费| 国语自产精品视频在线第100页| 黄色视频,在线免费观看| 欧美一级毛片孕妇| 日韩欧美国产一区二区入口| 日本一区二区免费在线视频| 欧美激情久久久久久爽电影| or卡值多少钱| 中文字幕久久专区| 免费高清视频大片| 可以在线观看毛片的网站| 免费av毛片视频| 高清在线国产一区| 免费av毛片视频| 国内久久婷婷六月综合欲色啪| 啦啦啦免费观看视频1| 精品一区二区三区四区五区乱码| 免费看a级黄色片| aaaaa片日本免费| 久久久久亚洲av毛片大全| 亚洲成人久久爱视频| ponron亚洲| 亚洲欧美日韩高清专用| 亚洲精品中文字幕在线视频| 在线观看免费午夜福利视频| 大型黄色视频在线免费观看| 亚洲人成伊人成综合网2020| 亚洲自拍偷在线| 亚洲一区高清亚洲精品| 亚洲成人中文字幕在线播放| 国产激情久久老熟女| 午夜a级毛片| 亚洲avbb在线观看| 精品国产美女av久久久久小说| 国产久久久一区二区三区| 草草在线视频免费看| 色尼玛亚洲综合影院| 草草在线视频免费看| svipshipincom国产片| 久久精品夜夜夜夜夜久久蜜豆 | 亚洲精品一区av在线观看| 国产日本99.免费观看| 亚洲真实伦在线观看| 99精品久久久久人妻精品| 免费在线观看影片大全网站| 午夜老司机福利片| 免费看a级黄色片| 国产69精品久久久久777片 | 免费人成视频x8x8入口观看| 老熟妇乱子伦视频在线观看| 色精品久久人妻99蜜桃| 欧美三级亚洲精品| 操出白浆在线播放| 男女午夜视频在线观看| 又爽又黄无遮挡网站| 丁香欧美五月| 欧美国产日韩亚洲一区| 丰满人妻一区二区三区视频av | 在线观看免费午夜福利视频| 人人妻,人人澡人人爽秒播| 两人在一起打扑克的视频| 天天躁狠狠躁夜夜躁狠狠躁| 欧美日韩亚洲国产一区二区在线观看| 在线永久观看黄色视频| 亚洲在线自拍视频| 亚洲无线在线观看| 欧美绝顶高潮抽搐喷水| 亚洲国产欧美一区二区综合| 精品国产亚洲在线| 别揉我奶头~嗯~啊~动态视频| 国产高清视频在线播放一区| av片东京热男人的天堂| 国产又色又爽无遮挡免费看| 国产aⅴ精品一区二区三区波| 99热只有精品国产| 麻豆一二三区av精品| 欧美性猛交黑人性爽| 国产亚洲精品av在线| 久久中文字幕一级| 免费无遮挡裸体视频| 男女视频在线观看网站免费 | 狂野欧美白嫩少妇大欣赏| 午夜免费成人在线视频| 久久久国产成人精品二区| 丁香欧美五月| 正在播放国产对白刺激| 日本精品一区二区三区蜜桃| 窝窝影院91人妻| 国产一区二区三区在线臀色熟女| 天天一区二区日本电影三级| 老司机靠b影院| 久久久国产成人免费| 久久精品人妻少妇| 成年女人毛片免费观看观看9| 国产精品香港三级国产av潘金莲| 观看免费一级毛片| 色综合站精品国产| 欧美绝顶高潮抽搐喷水| 一边摸一边抽搐一进一小说| 两个人的视频大全免费| 国内精品久久久久久久电影| 亚洲成人中文字幕在线播放| 日韩欧美国产在线观看| www.999成人在线观看| 欧美久久黑人一区二区| 午夜免费观看网址| 很黄的视频免费| 最近最新免费中文字幕在线| 熟妇人妻久久中文字幕3abv| 日韩欧美一区二区三区在线观看| 中出人妻视频一区二区| 日本免费一区二区三区高清不卡| 一本一本综合久久| 国产精品久久久久久久电影 | 精品不卡国产一区二区三区| 人成视频在线观看免费观看| 日韩国内少妇激情av| 亚洲人成77777在线视频| 可以在线观看毛片的网站| 国产精品一区二区三区四区久久| 色尼玛亚洲综合影院| 男女下面进入的视频免费午夜| 国产精品久久电影中文字幕| 99久久精品热视频| 中文在线观看免费www的网站 | 黑人欧美特级aaaaaa片| 国内揄拍国产精品人妻在线| 亚洲欧美精品综合久久99| 一进一出好大好爽视频| 天堂av国产一区二区熟女人妻 | 国产1区2区3区精品| 黄色视频,在线免费观看| 亚洲午夜精品一区,二区,三区| 久久精品成人免费网站| 日本免费一区二区三区高清不卡| 中出人妻视频一区二区| 88av欧美| 高潮久久久久久久久久久不卡| 国产精品日韩av在线免费观看| 一本综合久久免费| 这个男人来自地球电影免费观看| 91字幕亚洲| 色综合欧美亚洲国产小说| 宅男免费午夜| 91在线观看av| 欧美成人一区二区免费高清观看 | 国产一区二区激情短视频| 无限看片的www在线观看| 午夜福利在线观看吧| 一进一出抽搐gif免费好疼| 国产69精品久久久久777片 | 给我免费播放毛片高清在线观看| 婷婷精品国产亚洲av在线| 老司机福利观看| 国产成年人精品一区二区| 欧美日本视频| 一边摸一边抽搐一进一小说| 此物有八面人人有两片| 国产久久久一区二区三区| 亚洲精品美女久久av网站| 国产免费男女视频| 青草久久国产| 日韩欧美精品v在线| 丰满人妻熟妇乱又伦精品不卡| 国产av一区二区精品久久| 女人被狂操c到高潮| 亚洲成人中文字幕在线播放| 日韩精品免费视频一区二区三区| 亚洲成人久久爱视频| 不卡一级毛片| 欧美成人一区二区免费高清观看 | 国产精品永久免费网站| 日本免费一区二区三区高清不卡| 99国产精品一区二区蜜桃av| 亚洲欧美日韩高清在线视频| 亚洲国产看品久久| 热99re8久久精品国产| 欧美日韩瑟瑟在线播放| 亚洲熟妇中文字幕五十中出| 国产区一区二久久| 在线观看www视频免费| 男女之事视频高清在线观看| 日本在线视频免费播放| 精品久久久久久,| 在线观看免费日韩欧美大片| 9191精品国产免费久久| 宅男免费午夜| 蜜桃久久精品国产亚洲av| 亚洲精品在线观看二区| 18禁黄网站禁片午夜丰满| 波多野结衣高清无吗| 桃色一区二区三区在线观看| 村上凉子中文字幕在线| 女警被强在线播放| 男人舔女人的私密视频| 欧美丝袜亚洲另类 | 欧美成人午夜精品| 嫁个100分男人电影在线观看| 成人国产综合亚洲| 91av网站免费观看| 国产精品国产高清国产av| 色播亚洲综合网| 在线国产一区二区在线| 90打野战视频偷拍视频| 真人做人爱边吃奶动态| 久久精品影院6| www.精华液| 一a级毛片在线观看| 一进一出抽搐动态| 999久久久国产精品视频| 禁无遮挡网站| 亚洲专区国产一区二区| 亚洲中文av在线| 人妻丰满熟妇av一区二区三区| 一级a爱片免费观看的视频| 精品一区二区三区四区五区乱码| 色哟哟哟哟哟哟| 国产精品久久视频播放| 欧美性长视频在线观看| 久久久久久大精品| 九色国产91popny在线| 国产在线观看jvid| 国产欧美日韩一区二区精品| ponron亚洲| 久久久久性生活片| 精华霜和精华液先用哪个| 国产精品永久免费网站| 88av欧美| 日日摸夜夜添夜夜添小说| 丝袜美腿诱惑在线| 欧美不卡视频在线免费观看 | 久久久久国产精品人妻aⅴ院| 91成年电影在线观看| 黑人操中国人逼视频| 国产成人精品无人区| 欧美日韩瑟瑟在线播放| 久久久久九九精品影院| 99久久久亚洲精品蜜臀av| 久久久久久亚洲精品国产蜜桃av| 国内揄拍国产精品人妻在线| 精品久久久久久久久久免费视频| АⅤ资源中文在线天堂| 免费在线观看日本一区| 天堂√8在线中文| 中文字幕高清在线视频| 亚洲性夜色夜夜综合| 成人永久免费在线观看视频| 90打野战视频偷拍视频| 一本大道久久a久久精品| 亚洲黑人精品在线| 日韩中文字幕欧美一区二区| 久99久视频精品免费| 亚洲精品一区av在线观看| 欧美另类亚洲清纯唯美| 88av欧美| 亚洲午夜理论影院| 黄频高清免费视频| 激情在线观看视频在线高清| 给我免费播放毛片高清在线观看| 黄色成人免费大全| 久久午夜亚洲精品久久| 日韩精品青青久久久久久| xxxwww97欧美| videosex国产| 欧美乱码精品一区二区三区| 18禁美女被吸乳视频| 国产精品一区二区精品视频观看| 这个男人来自地球电影免费观看| 亚洲一区二区三区色噜噜| 亚洲人成伊人成综合网2020| 丁香欧美五月| 久久久久精品国产欧美久久久| 欧美黄色片欧美黄色片| 免费一级毛片在线播放高清视频| 国产欧美日韩精品亚洲av| 天堂av国产一区二区熟女人妻 | 麻豆成人午夜福利视频| 男人舔奶头视频| 国产日本99.免费观看| 日韩成人在线观看一区二区三区| 美女扒开内裤让男人捅视频| 国产成年人精品一区二区|