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

Online temperature estimation of Shell coal gasification process based on extended Kalman filter

2022-08-19 09:31:10KangchengWangJieZhangDexianHuang

Kangcheng Wang,Jie Zhang,Dexian Huang

1 Department of Automation,Beijing National Research Center for Information Science and Technology,Tsinghua University,Beijing 100084,China

2 School of Engineering,Merz Court,Newcastle University,Newcastle upon Tyne NE1 7RU,United Kingdom

Keywords:Shell coal gasification process Mechanistic modeling Temperature estimation Extended Kalman filter

ABSTRACT Obtaining the temperature inside the gasifier of a Shell coal gasification process (SCGP) in real-time is very important for safe process operation.However,this temperature cannot be measured directly due to the harsh operating condition.Estimating this temperature using the extended Kalman filter (EKF)based on a simplified mechanistic model is proposed in this paper.The gasifier is partitioned into three zones.The quench pipe and the transfer duct are seen as two additional zones.A simplified mechanistic model is developed in each zone and formulated as a state-space representation.The temperature in each zone is estimated by the EKF in real-time.The proposed method is applied to an industrial SCGP and the effectiveness of the estimated temperatures is verified by a process variable both qualitatively and quantitatively.The prediction capability of the simplified mechanistic model is validated.The effectiveness of the proposed method is further verified by comparing it to a Kalman filter-based single-zone temperature estimation method.

1.Introduction

With the depletion of oil resources,it is important to improve clean coal technology and use coal resources in a more efficient way,which also helps in reducing pollutant emissions [1].Coal gasification is one of the most important clean coal technologies,and the Shell coal gasification process(SCGP)is a popular coal gasification process.SCGP has advantages including large scales,high carbon conversion,and high cold gas efficiency over many other coal gasification processes[2].If SCGP can achieve stable operation for a long time,it is more profitable than most coal gasification processes.

The temperature inside the gasifier of SCGP has large impacts on gasification,syngas cooling,and slag formulation processes.An abnormal temperature will influence carbon conversion rate,cold gas efficiency,and syngas composition,and cause faults or even unscheduled plant shutdowns.Therefore,it is essential to obtain the temperature in real-time.However,considering the harsh environment in the enclosed gasifier,there is no temperature measurement device inside the gasifier[3].Consequently,this temperature cannot be obtained by operators.As a result,operators have difficulty responding quickly when the temperature changes abruptly,which leads to frequent faults and losses of economic benefits.

In recent decades,significant attention has been focused on developing mechanistic models of coal gasification processes and utilizing computer simulation to obtain the gasifier temperature.Wenet al.[5] analyzed the main chemical reactions in a Texaco downflow entrainment pilot plant gasifier.Based on material and energy balances,a mechanistic model was established and the temperature inside the gasifier was simulated.Watanabeet al.[6]established a model to describe three processes in an entrained flow coal gasifier:pyrolysis,char gasification,and gas-phase reactions.The simulated gas temperature distribution was presented.Monaghanet al.[7,31] developed a dynamic reduced order model to simulate entrained flow gasifiers.The simulation results of the model,including the syngas temperature,were validated by the simulation results of computational fluid dynamics and experimental records.Leeet al.[8] developed a detailed dynamic model for SCGP.The simulation data of the model,which contains the gas temperature,were validated by reference data.

Because the mechanistic models developed in these studies are rigorous,the computational costs of these models are huge[9].As a result,the gasifier temperature can hardly be obtained in realtime.To simplify the mechanistic models,the following temperature estimation models are developed.A multi-zone thermal calculation model for the pulverized coal boiler was established [10].A boiler was partitioned into several zones vertically with similar height.A heat balance model of each zone was built,and the gas temperature at the top of each zone was calculated by iteration.This model is efficient when high estimation accuracy is not required.However,this model is a static heat balance model,and the heat accumulation in the gasifier is not considered.Jiet al.[11] established simplified heat balance and mass balance equations of the syngas quenching and cooling processes of SCGP to calculate the temperature of syngas at the gasifier exit in real-time.However,the calculated temperature is not equal to the temperature inside the gasifier.Wanget al.[12]established simplified heat balance equations of SCGP.These equations were then formulated as a state-space representation.The state variable,the temperature inside the gasifier,was estimated by the Kalman filter (KF) algorithm.Nevertheless,this study assumed the temperatures in different regions of the gasifier are the same,while these temperatures are different in practice.

According to these studies,using simplified mechanism analysis to determine the model structure and utilizing parameter estimation methods to estimate the unknown parameters of the model is an effective way to obtain the gasifier temperature in real-time[13].Because the inaccuracies of the measurement devices can lead to unsatisfactory estimation results,an efficient noise filtration algorithm is needed to improve the estimation reliability [4].

The KF algorithm can estimate unknown variables by measurements containing noises [14,15].Wilsonet al.[16] utilized the KF algorithm to estimate the downstream pressure and feed coal quality of the ALSTOM gasifier.However,the gasifier temperature of the ALSTOM gasifier can be measured directly and served as an observation variable because it is much lower than that of the SCGP.More importantly,the system model developed in this work is linear.Considering SCGP is a complex system with strong nonlinearity,a nonlinear model together with a nonlinear temperature estimation algorithm are required to accomplish the estimation of gasifier temperature.

The extended Kalman filter (EKF) algorithm is an extension of the KF algorithm for state estimation of nonlinear systems [17].The EKF is very suitable for online estimation due to its computational efficiency.Recently,the EKF algorithm has also been applied to thermal systems whose time dynamics are modeled using heat flows.Mironovaet al.[18] utilized the EKF algorithm to estimate the temperatures of a Peltier cell and listed several advantages of estimating the temperatures by the EKF algorithm over measuring the temperatures by sensors.Valipouret al.[19] used the KF and EKF algorithms to estimate unknown states in entrained flow gasification systems based on a dynamic reduced order model.However,online estimation is still unachievable due to the relatively high complexity of the model.

For these reasons,a gasifier temperature estimation method is proposed in the present study.The gasifier of SCGP is divided into three zones.The quench pipe and the transfer duct are seen as two additional zones.A simplified mechanistic model is developed for each zone based on heat conservation.The mechanistic models are then formulated as a state-space representation whose state variables are the temperatures in individual zones.The EKF is adopted to implement the estimation of the temperature in each zone to obtain satisfactory estimation results with the existence of process noises and observation noises.Because the complexity of the mechanistic model is much lower than the existing mechanistic models of SCGP,the computational cost of the temperature estimation is very low,and the temperatures in different regions of the gasifier can be estimated in real-time.

The remainder of this article is arranged as follows:In Section 2,the main process units of SCGP are briefly introduced,and the main fault related to the temperature inside the gasifier is analyzed.In Sections 3 and 4,the gasifier is partitioned into several zones,and mechanistic models of individual zones are established to obtain a state-space representation.Then the EKF algorithm is applied to obtain the real-time estimation of temperatures in different zones,which can be validated by a process variable.In Section 5,the proposed temperature estimation method is applied to an industrial SCGP,and a detailed analysis of the temperature estimation results is carried out.The prediction capability of the simplified mechanistic model is validated,and a model comparison is conducted to further verify the effectiveness of the proposed method.At last,some conclusions are drawn in Section 6.

2.Process Description

The flow diagram of an SCGP is shown in Fig.1[20,21].It can be seen that the SCGP consists of several process units including feedstock preparation and feeding,gasification,syngas cooling,dry solids removal,and slag processing.

The feedstock of SCGP mainly contains coal,oxygen,and steam.The coal is dried and pulverized in advance,and then fed into the burner by high-pressure inert gases.The high-pressure oxygen from the air separation is preheated,mixed with the superheated steam,and fed in the gasifier through the burners.

In the gasifier,which is the main unit of the gasification process,feedstocks are burned intensely under the high-temperature and high-pressure condition.The main chemical reactions taking place in the gasifier together with their standard enthalpy of formation are as follows [22,23]:

The main components of the syngas produced by these chemical reactions are CO and H2.The syngas at the gasifier exit has a very high temperature (about 1550 °C).To prevent damage to devices in subsequent processes,the syngas will firstly be quenched in the quench pipe by the quench gas returning from the downstream.The quenched syngas will then be transferred through a transfer duct to the syngas cooler,where the syngas is further cool down.At the syngas cooler exit,the syngas has a relatively low temperature(about 250°C),which is bearable for downstream processes such as particulate removal.At last,the dry solids in the syngas is removed.Part of the clean cold syngas is recycled as the quench gas,and the remaining syngas is transferred to the downstream processes to be important raw materials [24].

The wall of the gasifier consists of several layers,and its structure is shown in Fig.2 [8].Among these layers,the membrane metal tube wall and the cooling water are designed to recover the heat,generate steam,and reduce the temperature of the gasifier wall.The other layers are employed to further reduce the wall temperature for reducing the heat loss to the ambient.Apart from the cooling water in the membrane metal tube wall,the cooling water in the evaporators of the quench pipe,the transfer duct,and the syngas cooler also absorbs a large amount of heat and generates steam.The steam is fed into the drum and be used in subsequent processes.

During the burning of coal,part of the ash in the coal forms slag.Some slag is in the solid state and sticks to the inner wall of the gasifier.Some slag is in the molten state.It flows down the inner wall and is discharged from the bottom of the furnace [25].Then it is cooled down and crushed.Both the solid-state and the molten-state slag can protect the wall from being destroyed by the high temperature inside the gasifier.

Many faults of SCGP are directly related to the temperature inside the gasifier.A typical example is the burnt-through of coal burners and burner covers.On one hand,the burners and burner covers tend to be burnt through if they work under a temperature higher than its normal range for a long period.This fault will lead to unscheduled plant shutdown and cause huge economic losses.On the other hand,if the temperature remains lower than its normal range for a long time,the fluidity of the slag on the inner wall of the gasifier will decrease,and a layer of slag will be formed at the entrance of the burner to block the flame sprayed by the burners.The blocked flame will burn the burner and burner covers in reverse.As a result,the burners and burner covers will also be burnt through.Therefore,the burners and burner covers are at risk either when the temperature is too high or too low.Therefore,obtaining the temperature in the gasifier in real-time is crucial for achieving the safe operation of SCGP.Since there is no temperature measurement device inside the gasifier,it is essential to build an online mathematical model to estimate the temperature.

3.Mechanism Analysis

To estimate the temperature inside the gasifier,the mechanism of an SCGP illustrated in Fig.1 is analyzed in this section.For simplicity,we assume the syngas temperature in each zone is of the same,and it is seen as the average temperature of the zone.

The gasifier is evenly partitioned into three zones along the height direction.The quench pipe and the transfer duct are seen as the fourth and the fifth zone.The partition scheme is illustrated in Fig.3.Noting that this figure is not drawn to scale.

By analyzing the mechanism of each zone,the following heat balance model can be derived [26,27]:

Fig.1.Flow diagram of SCGP.

Fig.2.Schematic of the structure of the SCGP gasifier wall.

Fig.3.Partition scheme of the gasifier,the quench pipe,and the transfer duct.

In addition,theTg,i(t)in Eq.(1)is equal to the syngas temperature at the entrance of the syngas coolerwheni=5,that is,

4.Temperature Estimation

4.1.Implementing the EKF

Assumingtsdenotes the sampling time of SCGP,then Eq.(1)can be transformed as follows by utilizing the explicit Euler method:

When considering process noises,Eq.(6) can be rewritten as follows:

where Tg=denotes the vector of temperatures in different zones;denotes the vector of nonlinear functions,which are expressed as the right side of Eq.(6);u=denotes the vector of input variables;v=denotes the vector of process noises,which are assumed to be white noises with zero means and autocorrelationwhere Q (k)is a symmetric nonnegative definite matrix [29].

Eq.(8)can be seen as the state equation of a state-space representation,and Tgdenotes a vector of state variables.The Jacobian matrix of f(·)in Eq.(8) is given as:

where

and the other elements of F (k)are equal to 0.

When considering observation noise,the following equation can be derived by Eq.(5):

wherey(k+1)=e(k+1)denotes the observation noise,which is assumed to be a white noise with zero mean and autocorrelationwhereR(k)is a positive scalar.Note that

Eq.(11)can be seen as the observation equation of a state-space representation,anddenotes the observation variable.

At last,by applying the EKF algorithm given as Eq.(12)to Eq.(8)and (11),the temperatures of different zones can be estimated in real-time [30].

4.2.Validating the estimation results

As stated in Section 1,the gasifier temperature of SCGP cannot be directly measured.As a result,the accuracy of the temperatures estimated by the EKF can hardly be calculated.Therefore,a measurable process variable is required to validate the effectiveness of the temperature estimation results.

When the reaction temperature increases,more heat will be absorbed by the cooling water in the membrane wall as well as the evaporators of the quench pipe and the transfer duct,and more steam will be generated in each evaporator.As a result,the mass flow rate of steam at the exit of the drumwill increase.Therefore,if the trend ofis consistent with the estimated temperatures,the estimation results can be considered as reasonable.Becausecan be measured online,and the time constant between the estimated temperatures andis short,is very suitable to be the validation variable.

Besides,the relationship between the estimated temperatures and thecan be described quantitatively by calculating the correlation coefficients between them.To reduce the interference of process noises,the estimated temperatures and theare firstly filtered by moving average filters,and then the correlation coefficient between the filtered values are calculated.

Note that the trend of the validation variable can also be used as additional information for operators to judge the gasifier temperature.

5.Application to an Industrial SCGP

In this section,an industrial SCGP is studied to verify the effectiveness of the proposed model.

5.1.Process description

The flow diagram of this process is shown in Fig.1.The measurements of all process variables are collected by a distributed control system (DCS).The sampling intervals of the process variables are 1 s,10 s,or 1 min,and they are unified as 1 min for modeling purpose.

In this process,the most commonly occurring fault is the burntthrough of burners and burner covers.Each fault leads to several days of unscheduled plant shutdown,causing huge economic losses.According to the fault analysis in Section 2,it is necessary to have an online estimation of the temperature inside the gasifier.The gasifier,the quench pipe,and the transfer duct are partitioned by the partition scheme in Section 3,and the temperatures in different zones are estimated in real-time by the EKF-based method in Section 4.

After removing outliers by the 3-sigma rule,a total of 2094 samples of process variables under normal operating conditions are selected for temperature estimation.The initial values of state variables are taken asbased on practical experience.The initial value of the covariance matrix in Eq.(12) is taken as P0(k)=100H5.The variances of processnoises and the observation noise are taken as Q (k)=100H5andR(k)=9,respectively.

The temperature estimation process is conducted on a Windows 10 computer with Intel Xeon CPU 3.3 GHz,8 GB RAM using MATLAB R2018b.The total time cost on the EKF estimation on all samples is around 0.4 s.

5.2.Results analysis

The estimated temperatures in the first three zones are shown in Fig.4(a)-(c).It can be seen that the estimated temperatures of these three zones have similar trends,which is reasonable because the temperatures in different zones of the gasifier are closely related.Besides,the values of the estimated temperatures decrease from zone 1 to zone 3,which is also reasonable because the burning percentage of fuel decreases as the height increases.In addition,the ranges of the three estimated temperatures decrease from zone 1 to zone 3,which is reasonable according to practical experience.

The estimated temperatures in the quench pipe and the transfer duct are shown in Fig.4(d)-(e).It can be seen that the trend of the estimated temperature in zone 4 is similar to that of the first three zones.This is because the ranges of the temperature and the mass flow rate of the quench gas are much lower than that of the syngas at the gasifier exit.As a result,the trend of the syngas temperature after quenching is similar to that of the syngas in the gasifier.The transfer duct is adjacent to the quench pipe,and therefore the trend of the syngas temperature in zone 5 is similar to that of zone 4.

It can also be seen that the variations in zones 1-3 are much severe than the variations in zones 4-5.This is reasonable because the syngas temperature becomes much more stable after quenching.Besides,the ranges of the estimated temperatures in zones 4-5 are much lower than that of the first three zones,which is also reasonable because the range of the syngas temperature decreases obviously after quenching.Moreover,the range of the estimated temperature in the transfer duct is lower than that of the quench pipe,which is consistent with practical experience.

To verify the effectiveness of the estimation results,the validation variable,i.e.,the mass flow rate of the steam at the drum exit,is illustrated in Fig.4(f).It can be seen that the estimated temperatures in all zones and the validation variable have similar trends,which is consistent with the mechanism analysis in Section 4.2.Besides,the variations in the validation variable are consistent with the variations in the estimated temperatures in zones 1-3,and higher than the variations in the estimated temperatures in zones 4-5.This is reasonable because the steam generated in the membrane wall of the gasifier is more than the steam generated in the evaporators of the quench pipe and the transfer duct.As a result,the temperatures in zones 1-3 have larger impact on the validation variable than zones 4-5.

The correlation coefficients between the estimated temperatures and the validation variable are calculated and listed in Table 1.It can be seen that the validation variable has very high correlation coefficients with the estimated temperatures in zones 1-3,which is reasonable because the validation variable is mainly influenced by the steam generated in the membrane wall.Besides,the validation variable also has relatively high correlation coefficients with the estimated temperatures in zones 4-5.Based on the above analysis,the estimated temperatures can be considered as being reasonable.

Table 1 Correlation coefficients between EKF estimated temperatures and the validation variable

Table 2 Correlation coefficients between predicted temperatures and the validation variable

In order to further verify the proposed temperature estimation method,a different time period of process operation data is collected for temperature estimation,and the results are shown in Fig.5.It can be seen that the results are similar to Fig.4.This means the proposed method is still effective on a different set of samples.

Fig.4. Estimated temperatures and the validation variable:(a)estimated temperature in zone 1;(b)estimated temperature in zone 2;(c)estimated temperature in zone 3;(d) estimated temperature in zone 4;(e) estimated temperature in zone 5;(f) mass flow rate of steam at the drum exit.

5.3.Model validation

In this section,the prediction capability of the simplified mechanistic model developed in Section 3 is validated.By artificially setting the Kalman gain to a zero matrix,the temperatures in different zones are predicted by the mechanistic model expressed as Eqn 6,and the prediction result is shown in Fig.6.By comparing Fig.6 to Fig.4,it can be seen that even without the correcting influence of the measurement innovation,the predicted temperatures in all zones are similar to the corresponding temperatures estimated by the EKF,and the trends of the predicted temperatures are consistent with that of the validation variable.Besides,the variations in the first three predicted temperatures are consistent with the variations in the validation variable.These results validate the prediction capability of the developed mechanistic model.

The correlation coefficients between the predicted temperatures and the validation variable are calculated and listed in Table 2.It can be seen that all correlation coefficients are close to the corresponding ones in Table 1.By comparing the last four correlation coefficients with the corresponding ones in Table 1,it can be seen that when the measurement innovation is incorporated in the prediction,the correlation coefficients in zones 2-5 are increased.This validates the effectiveness of the proposed EKF-based temperature estimation method.

Fig.5.Estimated temperatures and the validation variable on a different time period of process data:(a)estimated temperature in zone 1;(b)estimated temperature in zone 2;(c) estimated temperature in zone 3;(d) estimated temperature in zone 4;(e) estimated temperature in zone 5;(f) mass flow rate of steam at the drum exit.

Fig.6.Predicted temperatures and the validation variable:(a)predicted temperature in zone 1;(b)predicted temperature in zone 2;(c)predicted temperature in zone 3;(d)predicted temperature in zone 4;(e) predicted temperature in zone 5;(f) mass flow rate of steam at the drum exit.

Fig.7.(a) Predicted value of the observation variable;(b) measured value of the observation variable.

In addition,the predicted output of the observation equation is also calculated by the predicted temperatures through Eqn 5,and the calculation result is shown in Fig.7(a).The measured value of the observation variable,i.e.,the syngas temperature at the entrance of the syngas cooler,is shown in Fig.7(b).It can be seen that the predicted value of the observation variable is similar to the measured value.Besides,the predicted and measured values of the observation variable have similar trends.The correlation coefficient between the predicted and measured values of the observation variable is calculated as 0.7299.Moreover,the variations in the predicted value are consistent with the variations in the measured value.This further validates the prediction capability of the mechanistic model established in this study.

5.4.Model comparison

In this section,the proposed temperature estimation method is compared with the KF-based single-zone temperature estimation method proposed in Ref.[12].The gasifier temperature estimated by the latter method using the same process data is illustrated in Fig.8(a).The mass flow rate of steam at the drum exit shown in Fig.8(b) is also selected as the validation variable to validate the estimation result.It can be seen that the estimated temperature and the validation variable do not have similar trends on some samples.The correlation coefficient between the estimated temperature and the validation variable is calculated as 0.5401,which is much lower than that of the first three zones in Table 1.

Based on the above analysis,the temperature estimated by the method proposed in Ref.[12]is less reliable than the method proposed in this study.The main reason is the mechanistic model developed in Ref.[12] is relatively rough,and it has a relatively large error with the real process.

From the above model comparison,it can be seen that the temperature estimation method proposed in this article has the following two advantages over the method proposed in Ref.[12]:

(1) The proposed method can estimate the temperatures in different zones of the gasifier.As a result,the temperature near the burners can be estimated more accurately,which is crucial for operators to take immediate actions when the temperature goes out of its normal range.This will greatly reduce the fault rates of burners and burner covers.

(2) The temperature estimation results of the proposed method are much more reasonable.This means the proposed method can provide operators with more reliable estimated gasifier temperatures and help operators to make more accurate adjustments to operating conditions.

6.Conclusions

Since the temperature inside the gasifier of SCGP cannot be measured directly,online estimation of this temperature using the EKF based on a simplified mechanistic model is investigated in this paper.The gasifier is partitioned into three zones.The quench pipe and the transfer duct are seen as two additional zones.The mechanistic models in individual zones are established based on heat conservation.A state-space representation is developed based on the mechanistic models and the EKF is utilized to estimate the temperature in each zone.The estimation results of an industrial SCGP are verified by a process variable both qualitatively and quantitatively.The prediction capability of the simplified mechanistic model is validated.The proposed temperature estimation method is compared with the KF-based single-zone temperature estimation method,and the comparison demonstrates the effectiveness of the proposed method.

If the estimated temperatures go out of their normal ranges,operators can tune operating parameters immediately to bring the temperatures back to normal.In this way,the fault rate of gasifier temperature-related faults,especially the burnt-through of burners and burner covers,can be decreased.Therefore,estimating the temperatures is of great significance for extending the failurefree time length and increasing the economic benefits of SCGP.

Declaration of Competing Interest

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

Acknowledgements

The authors appreciate the funding from the National Natural Science Foundation of China (61673236 and 61873142) and the Seventh Framework Programme of the European Union (P7-PEOPLE-2013-IRSES-612230).The authors also thank Eng.Yue Yang,Prof.Fan Yang,Prof.Yongheng Jiang,Prof.Chao Shang,and Dr.Wenxiang Lv for helpful discussions.Kangcheng Wang especially appreciates the help of Dr.Haowei Hu.

Nomenclature

Aarea,m2

cpspecific heat capacity at constant pressure,kg·(kg·K)-1

Hspecific enthalpy,kJ·kg-1

Qheating value,kJ·kg-1

qpercentage of heat loss,%

Ttemperature,°C

ttime,s

Vvolume,m3

Δ D-value

δ fuel combustion fraction

ε emissivity

ρ density,kg·m-3

σ stefan-Boltzmann constant,W·(m2·K4)-1

ψ thermal efficiency coefficient

Superscripts

′at inlet

′′at outlet

amb ambient

dw downward

up upward

Subscripts

cal calculation

cf coal feed

cool cooling water

cs cross section

dm drum

g gas

gsf gasifier

ithei-th zone of gasifier,i=1,···,5

iw inner wall of gasifier

of oxygen feed

ph physical heat loss due to ash and slag

qch quench gas

s sampling

sgc syngas cooler

st steam

stf steam feed

uc unburnt carbon

ug unburnt gas

主站蜘蛛池模板: 国产成人无码AV在线播放动漫 | 精品91视频| 国产精品99在线观看| 97视频在线精品国自产拍| 国产尤物jk自慰制服喷水| 国产精品亚洲天堂| 极品私人尤物在线精品首页| 一区二区影院| 在线免费不卡视频| 日本爱爱精品一区二区| 国产精品福利导航| 国产午夜福利片在线观看| 亚洲男人的天堂在线| 免费一级毛片在线观看| 亚洲aⅴ天堂| 久久综合亚洲鲁鲁九月天| 国产在线一区二区视频| 久久精品视频亚洲| 男女猛烈无遮挡午夜视频| 亚洲an第二区国产精品| 欧美一区二区三区欧美日韩亚洲 | 成人国产一区二区三区| 91精品国产自产在线观看| 久久人人爽人人爽人人片aV东京热 | 中国毛片网| 91po国产在线精品免费观看| 狠狠躁天天躁夜夜躁婷婷| 毛片久久网站小视频| 久久a毛片| 成人免费一级片| 精品自窥自偷在线看| 色播五月婷婷| 久久久久九九精品影院 | 在线观看无码a∨| 国产白丝av| 青青青视频91在线 | 国产成人无码Av在线播放无广告| AV无码一区二区三区四区| 亚洲性视频网站| 国产第四页| 亚洲一欧洲中文字幕在线| 日韩在线观看网站| 免费一级α片在线观看| 亚洲第一成网站| 国产人成乱码视频免费观看| 超碰91免费人妻| 欧美日韩亚洲国产主播第一区| 国产精品亚洲精品爽爽| 在线看片中文字幕| 四虎永久在线视频| 国产精品欧美在线观看| 九色综合视频网| 久久成人免费| 91麻豆精品国产高清在线| 无码综合天天久久综合网| 九色国产在线| 亚洲天堂色色人体| 天天躁日日躁狠狠躁中文字幕| 国产美女免费| 精品久久国产综合精麻豆| 亚洲欧洲日产国产无码AV| 天堂网国产| 亚洲伊人电影| 亚洲码一区二区三区| 日本a∨在线观看| 欧美日韩导航| 亚洲娇小与黑人巨大交| 色网站在线视频| 亚洲天堂日韩在线| 91无码视频在线观看| 天天色天天操综合网| 日韩欧美91| аⅴ资源中文在线天堂| 亚洲自拍另类| 日本www色视频| 91精品视频播放| 亚洲成A人V欧美综合天堂| 日本欧美午夜| 久草视频中文| 免费大黄网站在线观看| 中文成人在线视频| 1024国产在线|