999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

Optimization for ASP flooding based on adaptive rationalized Haar function approximation☆

2018-09-28 07:08:52YuleiGeShurongLiXiaodongZhang

Yulei Ge ,Shurong Li,2,*,Xiaodong Zhang

1 College of Information and Control Engineering,China University of Petroleum(East China),Qingdao 266580,China

2 Automation School,Beijing University of Posts and Telecommunications,Beijing 100876,China

Keywords:Alkali-surfactant-polymer flooding Optimization Enhanced oil recovery Mathematical modeling Rationalized Haar function approximation Adaptive strategy

ABSTRACT This paper presents an adaptive rationalized Haarfunction approximation method to obtain the optimal injection strategy for alkali-surfactant-polymer(ASP) flooding.In this process,the non-uniform control vector parameterization is introduced to convertoriginal problem into a multistage optimization problem,in which a new normalized time variable is adopted on the combination of the subinterval length.Then the rationalized Haar function approximation method,in which an auxiliary function is introduced to dispose path constraints,is used to transform the multistage problem into a nonlinear programming.Furthermore,an adaptive strategy proposed on the basis of errors is adopted to regulate the order of Haar function vectors.Finally,the nonlinear programming for ASP flooding is solved by sequential quadratic programming.To illustrate the performance of proposed method,the experimental comparison method and control vector parameterization(CVP)method are introduced to optimize the original problem directly.By contrastive analysis of results,the accuracy and efficiency of proposed method are confirmed.

1.Introduction

Alkali-surfactant-polymer(ASP) flooding is an important tertiary oil recovery technology.The basic idea is:utilize three displacing agents(alkali,surfactant and polymer),whose synergistic effects are important to enhance oil recovery and to change the physicochemical property[1].Since the price of displacing agents is expensive,how to determine the injection strategy to maximize the profit as much as possible is always a challenging problem.

In current industrial application,the index comparison method,which is also called the experimental comparison method,is usually adopted to select the injection strategy,in which the best one is chosen among some given feasible strategies by simulating on the numerical simulation software according to a defined index[2].However,it is too dependent on the experience of experts.The optimization result may not be optimal.While the optimal control technology can search for the optimal control strategy from all feasible solutions with considering the process dynamic characteristics.It can realize the simultaneous optimization of all variables.It is suitable for the optimization of ASP flooding.Some scholars have studied on the optimal control for oil exploitation.Lei[3]presented a mixed-integer iterative dynamic programming(IDP)to optimize the polymer flooding and had gotten good result.Ramirez and Fathi[4]optimized the injection strategy for enhanced oil recovery of surfactant flooding with optimal control theory,in which the necessary condition of optimal control was deduced on the basis of maximum principle before getting solved by gradient method.Furthermore,this method was applied to carbon dioxide flooding,nitrogen flooding and binary system flooding[5].Most of the researches are about water flooding and single chemical flooding.The optimal control for ASP flooding needs to be further studied urgently.

Numerical methods are usually adopted in optimization,which can be divided into two classes:direct methods and indirect ones.Direct methods include the control vector parameterization(CVP)[6],the direct multiple shooting[7],the full discretization[8],etc.What differs with these methods is how to select the kind of variables to be discretized and how to approximate the systemequations.However,direct methods may bring about unsatisfactory results especially when the control variables are discontinuous,such as with switching points and singular arcs[9].For these reasons,the reasonable discretization methods for variables and approximation methods for system equations are quite important.On account of the above,an adaptive CVP method was proposed in reference[10],in which the problem was discretized adaptively over time spans.The discretization was sequentially refined based on a wavelet-based analysis of the optimal solution which was obtained in the previous optimization step.As to the approximation methods for system equations,Haar function approximation method,which is a major kind of direct methods,has been widely used in recent years[11].Li[12]presented a numerical method based on quasilinearization and rationalized Haar function to solve nonlinear constrained optimal control problems.Oru?[13]investigated numerical solutions of the regularized long wave equation by using Haar wavelet,combined with finite difference method.Swaidan[14]proposed a numerical method for solving nonlinear optimal control problem with state and control inequality constraints.After being processed by the Haar wavelet operational matrix,the method was applied to solve optimal control of multi-item inventory model.Haar function,which only consists of three values:+1,?1 and 0,has simple structure and is easy to use.Although many researchers have studied on Haar function methods,the order of operational matrix in most of them cannot regulate adaptively.In calculation,if the order is too low,the accuracy cannot meet our demand.But if it is too high,this will increase computational burden.It is important to find a suitable order.

In view of above presentations,an optimization method for solving the optimal injection strategy for ASP flooding is developed in this paper.The original optimization problem is transformed into a multistage problem by non-uniform CVP firstly.Then the adaptive rationalized Haar function approximation is adopted to approximate the states and controls with errors regulating the order of operational matrix.After being transformed into a nonlinear programming(NLP),the problem is solved by sequential quadratic programming(SQP).Finally,the experimental comparison method(ECM)and CVP method are introduced as comparison to optimize the enhanced oil recovery of ASP flooding to illustrate the effectiveness of proposed method.

2.Problem formulations

2.1.Mechanism model description for ASP flooding

In the optimal control problem for ASP flooding,the control variables are the injection concentrations of three displacing agents(alkali,surfactant and polymer)at all injection wells,the states cover grid concentration and water saturation,and the performance index is the maximal NPV.This is a complex distributed parameters control problem aiming to get the optimal injection strategy to fulfill the maximal profit and enhanced oil recovery[15].Since the three-dimensional model for ASP flooding,which includes a series of divergence and cross terms,is too complex,only the one-dimensional model is considered in this paper.

Considering a long tube core with diameter d and length L,inject the displacing agents liquid with flow Q from one side at a constant speed,the core porosity is ? and the residual oil saturation is Sor.On the basis of the model assumptions in references[1,2],we add the below assumptions:

(1)The stratum satisfies heterogeneity,the whole reservoir is isothermal,and the adsorption process complies with the Langmuir isothermal adsorption equation;

(2)The oil displacement system is composed of alkali,surfactant and polymer.All displacing agents exist in water phase,while only surfactant exists in oil phase;

(3)The Darcy law is fit for the flow of oil phase and water phase;

(4)The phase equilibrium is set up instantly for all kinds of adsorptions which satisfy the generalized Fick's law;

(5)The fluid and rock are slightly compressible.Consider the effects of capillary force and gravity,and ignore the quality variation caused by the addition of displacing agents.Consider the permeability changes and the inaccessible pore volume of polymer.

On the basis of the oil/water seepage continuity equations and adsorption diffusion equations of displacing agents,combining mass balance conditions we can obtain below model with interaction of alkali,surfactant and polymer fully considered.

The seepage continuity equation is

The adsorption diffusion equations of surfactant and polymer are

The adsorption diffusion equation of alkali is

where a denotes the alkali,s denotes the surfactant,p denotes the polymer,Swis the water saturation,A is the core cross section area,fwis the moisture content,vwdenotes the seepage speed of water phase,Ca,Cs,Cpand Da,Ds,Dpdenote the concentration and diffusion coefficient of alkali,surfactant and polymer respectively.ρrdenotes the core density,Γa,Γs,and Γpare the adsorbing capacity of core for different displacing agents,and Rais the alkali consumption.

The initial conditions are

The boundary conditions are

where uΘdenotes the injection concentration of displacing agents,which is the control variable.

In application,the slug injection strategy is usually adopted.Suppose that there are P slugs.

where tidenotes the time node,the length of every slug is Ti=ti? ti?1,uΘidenotes the injection concentration of displacing agent Θ in slug i.Fig.1 shows a three-slug injection scheme.

Furthermore,the dosage limit of displacing agents is

where MΘpdenotes the maximum usage of displacing agents Θ.

The injection concentration constraints and slug size constraints are

The maximum economic profit is chosen as the index,the specific description is

where χ denotes the discountrate,qin,qoutdenote the volume flow rate of injection and production.

Fig.1.Three-slug injection scheme.

The other physicochemical algebraic equations can be found in Appendix A.

2.2.Optimization problem formulations

To make the problem more concise,we use w ∈ Rnw(CΘ,Sw,ΓΘ)to denote all states and u∈Rnu(the injection concentration of all displacing agents uΘ)to denote all control variables.Apply the finite difference method[19]to discretize the space variable,and keep the time variable same.Thus the PDEs can be transformed into ODEs with respect to time t.Then the optimization for ASP flooding can be reformulated as below problem.

where f(g)∈Rnfdenotes the system process function.The path constraints and terminal constraints are expressed as c(·)∈ Rncand φ(·) ∈ Rnφ,respectively.uLB,uUBdenote the lower bound and upper bound of control.

3.Adaptive rationalized Haar function approximation

3.1.Multistage problem formulation

Since the three-slug strategy is applied in oil exploitation,the control variable in Eqs.(10a),(10b),(10c),(10d),(10e)is discontinuous.To cope with the discontinuity points of control variables preferably,the non-uniform CVP is adopted to convert the original problem into a multistage problem(MSP)[16].The main idea of this process is to approximate the control with piecewise function.Discretize the whole time domain[t0,tf]into K uneven time periods.The specific time control nodes are

Approximate the continuous control variable with piecewise function in every period,that is

Furthermore,a new time variable τ is defined to deal with the time variable,

Then every original subinterval t∈[tk?1,tk]is transformed into τ∈[0,1].

Take the derivation of Eq.(13)with respect to τ,then

Define w(k)(τ),u(k)(τ)as the state and control on k th period,respectively.Then Eqs.(10a),(10b),(10c),(10d),(10e)can be reformulated as below.

Eq.(15e)is the linkage constraint which is used to keep the continuity of state variables at the interface of all subintervals.Eq.(10e)can be merged into Eq.(10d)and is transformed into Eq.(15d).The MSP in Eqs.(15a),(15b),(15c),(15d),(15e)consists of K sub-problems corresponding to K subintervals.The final control is obtained by solving K optimal control problems.

3.2.Rationalized Haar function approximation

Once the MSP as shown in Eqs.(15a),(15b),(15c),(15d),(15e)is obtained,the next step is to transfer it into a discrete NLP,which can be solved numerically by well-developed algorithms.In this section,this transformation is achieved by rationalized Haar function approximation[13].On the basis of the presentation about Haar function approximation in reference[12],below deductions are made.

By applying the rationalized Haar function to each state and control in subinterval k,the state and control can be approximated as

where In,Imdenote the n×n,m×m positive definite matrix.?denotes the Kronecker product.

The initial condition for subinterval k is

Here d is a nr×1 vector.

Integrate Eq.(15b)from tk?1to tk,below equation can be obtained

The system equations in Eq.(15b)can be written as

variables are used for both w(k)(1)and w(k+1)(0),constraints of Eq.(15e)can be removes and constraints of Eq.(22)can be written as

where 1≤k≤K?1.

The inequality constrains in Eq.(15d)can be transformed into equality constraints by introducing ancillary function z(τ).

z(τ)can be approximated rationalized by Haar function as z(τ)=τ)Z.Then Eq.(15d)can be reformulated as

Then another optimization sub-problem can be got,

Similarly,Eq.(15c)can be reformulated as

The performance index function of Eq.(15a)can be approximated with the similar method as Eq.(22).

Finally,the NLP that arises from rationalized Haar function approximation is used to maximize the performance index function of Eq.(29),subjecting to algebraic constraints of Eq.(23a),(23b),(25a),(25b)and(27).

Optimal solution of the resulting NLP fulfills the KKT conditions of optimality based on the augmented performance index which is defined as Eq.(30).λ,μ,and v are the Lagrange multipliers associated with process descriptions of Eq.(23a),(23b),path constraints of Eq.(25a),(25b),and terminal constraints of Eq.(27)in subinterval k,respectively.The KKT optimality conditions can be acquired by differentiating Eq.(30)with respect to every variable and setting them to zero.

3.3.Adaptive strategy

In the process of rationalized Haar function approximation,the accuracy of results mainly depends on the order of operational matrix.The order has to be chosen high enough to meet the given accuracy,although this will arouse the burdensome calculation.To regulate the order adaptively,an adaptive strategy is proposed in this section to obtain a reasonable order,which can balance the accuracy and computational burden.

In the subinterval k,the estimation of all statescan be obtained by Eq.(21).Define the absolute error and relative error as follows.

where i=1,2,K,n.

Then the maximal relative error can be written as

Assume that the given accuracy tolerance is ε.For subinterval k,ifwe can conclude that the current order is reasonable.Ifwe can adopt the below equation to update the order[17].

Here,ceil is the operator that rounds to the next highest integer.

On the other hand,assume that rmin,rmaxare the lower bound and upper bound of the order of operational matrix for rationalized Haar function approximation in subinterval k.If the regulated ordersatisfies that≤rmax,thenis adopted to update the order.On the other hand,if>rmaxwhich means computational burden is too large,then the current subintervals should be further refined.Update the number of subintervals with K′.

4.Optimization procedures

In this section,the SQP[18]is adopted to solve the NLP that is described as Eqs.(23a),(23b),(25a),(25b),(27)and(29).The specific procedures are as follows.

(1)Initialize the parameters,rmax,rmin,the number of subintervals K,the initial state w(1)(0)=0,initial control U(1)(0)=0,the accuracy tolerance εJ,ε,and the current iteration l=1;

(2)Transcribe original problem to NLP with current subintervals;

(3)Calculate the states and performance index Jlof each subinterval with current order for operational matrix in step l according to Eqs.(21)and(30);

(6)Assess the refinement of subintervals:If≤Nmax,go to step 7),if not,update the subintervals as Eq.(35);

(7)Determine the search direction and search step length,solve the NLP,and calculate the optimal controland index Jl+1;

If|Jl+1?Jl|<εJ, finish the computing process and output the optimal control Ul+1and index Jl+1.If not,go to step 2)and let l=1,carry on the procedure until|Jl+1?Jl|satisfies the given accuracy.

5.Optimization result for ASP flooding

In this section,the optimization problem of ASP flooding in Section 2.2 is solved by the proposed method in this paper.The three slug injection is adopted to calculate.There are twelve optimization variables which are uΘ1,uΘ2,uΘ3,TΘ,where TΘis the slug size for all displacing agents.

Some essential data for solving the optimization model is shown in Table 1.

Table 1 Partial data for solving ASP flooding model

Before optimization,we use the experimental comparison method(ECM)to select a relatively optimal solution artificially.Give lots of different control inputs,and compute the performance index.By the comparison of many times experiment results,we select the best injection concentration and slug size with the maximum profit.The injection strategy is

Then transform the ASP flooding problem into a NLP with proposed method which is presented in Section 3.The parameters are set as follows.

To illustrate the accuracy of proposed method,non-uniform CVP method is introduced to solve the original problem in Section 2.2.After being optimized,the optimal injection strategy with CVP is

The dimensionless profit is J?=0.26385.The running time is 126.36 s.The optimal injection strategy with proposed method is

The dimensionless profit is J?=0.26459.The running time is 79.76 s.The results are shown in Figs.2–5.

Fig.3.Comparison results of injection concentration for surfactant.

Fig.4.Comparison results of injection concentration for polymer.

Fig.5.Comparison results of performance index.

According to these figures,the index of proposed method is 0.01079 more than that of ECM,and is very close to that of CVP.So the result of proposed method is good and is consistent with the theoretical solution,basically.Furthermore,the running time of proposed method is less than that of CVP method.Since the usage of adaptive strategy,the order of operational matrix which can balance the error and computational time is determined,and the discontinuity points are processed more intensively.Since the proposed method has good accuracy and operating efficiency,it is suitable for the optimization of ASP flooding.

6.Conclusions

In this paper,a method based on adaptive rationalized Haar function approximation has been proposed to solve optimal injection strategy for ASP flooding whose control variables are discontinuous.Firstly,the original problem is converted into a multistage problem with non-uniform CVP method.Then the rationalized Haar function method was used to approximate the states,controls and system process,and to transform the optimization for ASP flooding into a NLP.To regulate the discontinuity points and balance the accuracy and efficiency,an adaptive strategy was developed based on errors.At last,the problem was optimized with SQP.Numerical simulations were carried out to show the effectiveness of the proposed method compared with CVP and ECM.The method in this paper can deal with discontinuity points well,it can be used to solve one kind of optimization problem with discontinuous control variables.

Nomenclature

Subscripts

References[1]H.B.Li,The New Progress and Field Test Research for ASP Flooding,Science Press,Beijing,2007 115–149.

[2]S.R.Li,X.D.Zhang,Optimal Control of Polymer Flooding for Enhanced Oil Recovery,China University of Petroleum Press,Dongying,2013 45–99.

[3]Y.Lei,S.R.Li,X.D.Zhang,Q.Zhang,L.L.Guo,Optimal control of polymer flooding based on mixed-integer iterative dynamic programming,Int.J.Control.84(11)(2011)1903–1914.

[4]W.F.Ramirez,Z.J.Fathi,C.Cagnol,Optimal injection policies for enhanced oil recovery:Part l–theory and computational strategies,Soc.Pet.Eng.J.24(3)(1984)328–332.

[5]W.F.Ramirez,Application of Optimal Control to Enhanced Oil Recovery,Elsevier,Amsterdam,1987.

[6]V.Vassiliadis,R.Sargent,C.Pantelides,Solution of a class of multistage dynamic optimization problems:1.Problems without path constraints,Ind.Eng.Chem.Res.33(1994)2111–2122.

[7]U.Ali,Y.Wardi,Multiple shooting technique for optimal control problems with application to power aware networks,IFAC 48(27)(2015)286–290.

[8]L.T.Biegler,Solution of dynamic optimization problems by successive quadratic programming and orthogonal collocation,Comput.Chem.Eng.8(1984)243–247.

[9]B.Srinivasan,S.Palanki,D.Bonvin,Dynamic optimization of batch processes I.Characterization of the nominal solution,Comput.Chem.Eng.27(2003)1–26.

[10]M.Schlegel,K.Stockmann,T.Binder,W.Marquardt,Dynamic optimization using adaptive control vector parametrization,Comput.Chem.Eng.29(2005)1731–1751.

[11]M.Mansoori,A.Nazemi,Solving in finite-horizon optimal control problems of the time-delayed systems by Haar wavelet collocation method,Comput.Appl.Math.2014(2014)1–21.

[12]Z.Y.Han,S.R.Li,Quasilinearization and rationalized Haar approach for direct solution of optimal control problems,Corpoica Cienc.Tecnol.Agrop.14(9)(2013)511–516.

[13]?.Oru?,F.Bulut,A.Esen,Numerical solutions of regularized long wave equation by Haar wavelet method,Mediterr.J.Math.2016(2016)1–19.

[14]W.Swaidan,A.Hussin,Haar wavelet operational matrix method for solving constrained nonlinear quadratic optimal control problem,AIP Conf.Proc.1682(1)(2015)69–72.

[15]Y.L.Ge,S.R.Li,P.Chang,R.L.Zang,Y.Lei,Optimal control for an alkali/surfactant/polymer flooding system,The 35th Chinese Control Conference,2016 2016,pp.2631–2636.

[16]P.Wang,C.H.Yang,Z.H.Yuan,The combination of adaptive pseudospectral method and structure detection procedure for solving dynamic optimization problems with discontinuous control profiles,Ind.Eng.Chem.Res.53(17)(2014)7066–7078.

[17]C.L.Darby,W.W.Hager,A.V.Rao,Direct trajectory optimization using a variable low-order adaptive pseudospectral method,J.Spacecr.Rocket.48(2011)433–445.

[18]A.F.Izmailov,A.S.K.M.V.Solodov,Some composite-step constrained optimization methods interpreted via the perturbed sequential quadratic programming framework,Optim.Methods Softw.30(3)(2015)461–477.

[19]M.Marchewka,Finite-difference method applied for eight-band kp model for Hg1?xCdxTe/HgTe quantum well,Int.J.Mod.Phys.B 31(20)(2017),1750137.[20]C.Z.Yang,et al.,Enhanced Oil Recovery for Chemical Flooding,Petroleum Industry Press,Beijing(China),2007.

[21]Y.L.Ge,S.R.Li,P.Chang,S.L.Lu,Y.Lei,Optimization of ASP flooding based on dynamic scale IDP with mixed-integer,Appl.Math.Model.2017(2017)1–16.

[22]Y.L.Ge,S.R.Li,K.X.Qu,A novel empirical equation for relative permeability in low permeability reservoirs,Chin.J.Chem.Eng.22(11-12)(2014)1274–1278.

Appendix A

The relevant physicochemical algebraic equations are shown as below:

A1.Alkali consumption

Many kinds of conditions can influence the alkali consumption,while only the adsorption consumption,acidoid consumption and heavy metal ion consumption are considered in this paper to describe the mechanism well and simply.The alkali consumption is given through the chemical reaction term ri,which is defined as follows:

where riis the alkali consumption per unit volume of different effects,the specific forms are as follows:[20]

(1)The rapid alkali consumption r1caused by the ion exchange between Na+in solution and H+on rock surface.This effect can be expressed as the Langmuir isothermal adsorption equation:

where r1is the maximal alkali consumption of this effect,a1is the coefficient,they are set by experimental data.

(2)The alkali consumption r2caused by acidoid is

The curve of r2under different acidity and alkalinity can be obtained through experiment.

(3)The alkali consumption r3caused by Ca2+,Mg2+acidoid in oil is

where KCa2+,KMg2+denote the solubility product.

A2.Surface tension

Surface tension varies with the concentration of alkali,surfactant and polymer.The form is[20]

A3.Capillary pressure

The capillary pressure of compound system can be described as the function of oil/water capillary pressure and surface tension,that is[20]

where CPC,NPCdenote the constant,Snis the wet phase saturation.

A4.Surfactant adsorption

Considering that the relation of polymer adsorption is reversible with salinity,and is irreversible with concentration,the adsorption can be described as follows according to Langmuir isothermal adsorption law[21].

A6.Relative permeability

The relative permeability of oil phase kroand water phase krware usually given by interpolation method in the process of model solution.Here the formula proposed in reference[22]is adopted.

where A,B,C,and D denote the coefficients which are identified by the data of kroand krw.

Since the addition of displacing agent,the permeability reduction factor has to be considered.

where qθ,Rkdenote the adsorption capacity of displacing agent and permeability reduction factor of water phase under different salinity,respectively.

A7.Liquid viscosity

For the displacing liquid of ASP flooding,the relation between viscosity μ and shear rate˙γ can be expressed by Meter equation as[21]

A8.Residual saturation

The residual saturation is relevant to the capillary number Nc,and can be obtained through the experiments[21].

where Srjdenotes the residual saturation of phase j.

主站蜘蛛池模板: 72种姿势欧美久久久大黄蕉| 国产精品人人做人人爽人人添| 亚洲第一极品精品无码| 国产欧美日韩一区二区视频在线| 亚洲精品欧美日韩在线| a毛片免费观看| 国产白浆一区二区三区视频在线| 五月婷婷综合网| 丰满人妻久久中文字幕| 国产高潮视频在线观看| 国产一区二区人大臿蕉香蕉| 一区二区自拍| 亚洲中文字幕国产av| 国产激情无码一区二区APP| 伊在人亚洲香蕉精品播放| 九九久久精品国产av片囯产区| 婷婷综合缴情亚洲五月伊| 亚洲一区色| 91精品国产综合久久香蕉922| 色妞永久免费视频| 国产福利微拍精品一区二区| 日韩无码真实干出血视频| 国产精品尹人在线观看| 日韩在线第三页| 国产精品久线在线观看| 欧美国产综合色视频| 亚洲成人在线免费| 亚洲香蕉伊综合在人在线| 亚洲色图综合在线| 国产福利免费视频| 欧美亚洲中文精品三区| 成人在线亚洲| 国产精品一区二区国产主播| 久久久久久尹人网香蕉 | 国产成人av大片在线播放| 香蕉视频在线精品| 天天综合网亚洲网站| 国产成人免费观看在线视频| 国产日本一区二区三区| 欧美区日韩区| 亚洲无线视频| 亚洲精品视频免费看| 日本午夜三级| 欧美在线伊人| 国产成人一区在线播放| av无码久久精品| 亚洲色无码专线精品观看| 91丝袜在线观看| 日韩久久精品无码aV| 日韩无码真实干出血视频| 91精品国产福利| 亚洲人成色在线观看| 丰满人妻被猛烈进入无码| 国内精品免费| 日本欧美在线观看| 四虎影视库国产精品一区| 精品福利视频导航| 亚洲人成影院在线观看| 国内毛片视频| 四虎永久在线视频| 午夜a级毛片| 欧美精品伊人久久| 亚洲欧美日韩久久精品| 制服丝袜国产精品| 亚洲欧美h| av尤物免费在线观看| 国产午夜福利片在线观看| 精品久久国产综合精麻豆| 91麻豆久久久| 青青青草国产| 毛片免费高清免费| 欧美不卡二区| 伊人色天堂| 欧美日韩国产成人高清视频| 亚洲第一视频免费在线| 黄色在线网| 久久国产拍爱| 久久免费观看视频| 亚洲成A人V欧美综合| 伊人成人在线| A级全黄试看30分钟小视频| 日日碰狠狠添天天爽|