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

Study on gas–liquid flow characteristics in stirred tank with dual-impeller based on CFD-PBM coupled model

2021-12-08 13:31:16SongsongWangQiuxiangBuDeyuLuanYingZhangLongbinLiZhaoruiWangWenhaoShi
Chinese Journal of Chemical Engineering 2021年10期

Songsong Wang,Qiuxiang Bu,Deyu Luan,Ying Zhang,Longbin Li,Zhaorui Wang,Wenhao Shi

College of Electromechanical Engineering,Qingdao University of Science and Technology,Qingdao 266061,China

Keywords:Gas–liquid two-phase flow Gas holdup Gas–liquid interfacial area Bubble Coupling algorithm of CFD-PBM

ABSTRACT Study on gas–liquid flow in stirred tank with two combinations of dual-impeller(six-bent-bladed turbine(6BT)+six-inclined-blade down-pumping turbine (6ITD),the six-bent-bladed turbine (6BT)+six-inclinedblade up-pumping turbine(6ITU))was conducted using computational fluid dynamics(CFD)and population balance model (PBM) (CFD-PBM) coupled model.The local bubble size was captured by particle image velocimetry(PIV)measurement.The gas holdup,bubble size distribution and gas–liquid interfacial area were explored at different conditions through numerical simulation.The results showed that the 4 mm bubbles accounted for the largest proportion of 33% at the gas flow rates Q=0.76 m3.h-1 and 22% at Q=1.52 m3.h-1 for combined impeller of 6BT+6ITU,while the bubbles of 4.7 mm and 5.5 mm were the largest proportion for 6BT+6ITD combination,i.e.25% at Q=0.76 m3.h-1 and 22% at Q=1.52 m3.h-1,respectively,which indicated that 6BT+6ITU could reduce bubble size effectively and promote gas dispersion.In addition,the gas holdup around impellers was increased obviously with the speed compared with gas flow rate.So it was concluded that 6ITU impeller could be more conductive to the bubble dispersion with more uniform bubble size,which embodied the advantages of 6BT+6ITU combination in gas–liquid mixing.

1.Introduction

The gas–liquid stirred tanks are widely used in industrial production owing to its high heat and mass transfer efficiency,large interphase area and simple operation [1–3].Many studies have been devoted to the gas–liquid flow characteristics [4–8],overall gas holdup distribution [9,10],local gas holdup [11–13],power change with gas flow rate [14,15],and distributions of bubbles sizes [16–18].In recent years,computational fluid dynamics(CFD) simulation has been applied to study the gas–liquid twophase flow characteristics in a stirred tank by a growing number of researchers.Min et al.[19]studied the gas void fraction distribution in stirred tank with blade disk turbine and two up-pumping maxflow impellers (BT-6+2MFU) combination impeller through a population balance model combined with a multiple size group(MUSIG)model.The results showed the MUSIG model could accurately predict the local gas holdup above the level of the top impeller,which was in good agreement with the experimental results.Bao et al.[20] investigated the effects of the impeller diameter on gas dispersion with parabolic blade disk turbine and two down-pumping hydrofoil propellers (PDT+2CBY) combination by coupling algorithm of CFD and PBM (population balance model)model.Their results showed that the gas holdup increased significantly with impeller diameter at lower superficial gas velocity.Yang et al.[21] analyzed the effects of gas flow rate on the flow field structure,gas holdup and power consumption using dualimpeller composed of standard Rushton and dislocated-blade Rushton,respectively,by CFD simulations.They found that dualimpeller of dislocated-blade Rushton had a certain superiority at same operating conditions through the comparison with the experimental results.Ranade et al.[22]investigated the influence of gas flow rate on the structure of trailing vortices for a Rushton turbine based on computational snapshot approach.The results indicated that the computational model appeared to capture the key characteristics of the gas–liquid flows in stirred tanks,including the gas cavern behind impeller blades.Aubin et al.[23]explored the effects of aeration on mean liquid velocities.They found that the presence of gas only slightly altered the liquid flow patterns produced by both the down and up-pumping impellers and caused a general decrease in the mean liquid velocities.Zhu et al.[24] studied the gas–liquid mixing in a bioreactor equipped with an up-or a down-pumping impeller under unaerated and aerated conditions.The results showed that the up-pumping mode offered the apparent advantages over the down flow.Zhang et al.[25,26]investigated the effects of coalescence and breakup models on gas holdup and local bubble size distribution in a bubble column and the effects of four parameters (rv,dmin,dmaxand dc) on local gas holdup in the Class Method for solving PBEs.They found that Luo’s coalescence and breakup model could predict the gas holdup well by modifying coalescence rates with a factor ranging from 0.4 to 0.8.They proposed an optimization scheme in the Class Methods(CM) model.Gu et al.[27] investigated the influence of impeller speed and superficial gas velocity on gas dispersion in stirred tanks with three combinations of multiple impellers based on the CFDPBM model.The results indicated that punched rigid-flexible combination impeller was more efficient in terms of gas dispersibility than rigid and rigid-flexible combination impeller.

It is well known that the multiple-impeller is usually applied in the mixing operation owing to the increasing height diameter ratio.The studies on gas–liquid dispersion in stirred tank with multiple-impeller have indicated that the connecting flow pattern can be formed between the impellers when a radial flow impeller is located at the bottom and an axial flow impeller is installed at the above,which can enhance the overall circulation of fluid and facilitate the dispersion of bubbles [28,29].From the literature review,it is clear that the mixing performance of impellers in gas–liquid two-phase flow needs to be investigated as little research has been devoted to the effects of impeller style on the gas–liquid dispersion by the numerical simulation.Therefore,the aim of this work is to study the gas–liquid mixing characteristics in stirred tank with two combinations of dual-impeller (six-bentbladed turbine (6BT) as the lower impeller,six-inclined-blade down-pumping turbine (6ITD) and six-inclined-blade uppumping turbine (6ITU) as the upper impeller,respectively,i.e.6BT+6ITDand 6BT+6ITU) comparatively based on the CFD-PBM coupled model.The distribution of gas holdup,bubble size and gas–liquid interfacial area were explored in detail.The results showed that combined impeller of 6BT+6ITUhad an obvious advantage on the gas–liquid mixing,which can provide reference for the optimization of dual-impeller and the application of gas–liquid two-phase mixing in industry.

2.Experimental

2.1.Experimental setup

The experimental setup of gas–liquid mixing is shown in Fig.1.The main devices include the lifting frame,stirred equipment,and electromagnetic speed regulating motor and data collection system.

Fig.1.The experimental setup of gas–liquid mixing.

The PIV experiment apparatus is shown in Fig.2,which mainly includes the encoder,laser transmitter,image acquisition system,synchronizer controlling system,computer,and image data processing system.The image acquisition system is composed of two cross-frame CCD cameras with a resolution of 2048 × 2048 pixels.The laser sheet thickness and laser pulse interval was adjusted to 1 mm and 220 μs,respectively.The sampling frequency was all set on 8 Hz.

The gas–liquid stirred equipment consists of a cylindrical tank(dia.Ts=420 mm) with a standard ellipsoidal head at the bottom of tank,as shown in Fig.3.Four equally spaced baffles (width:w=Ts/10,height:hb=400 mm)were fitted in the tank,The clearance from baffle to the tank wall is 2 mm.Tap water(density:ρl=998.2 kg.m-3,dynamic viscosity:μl=1 × 10-3Pa.s) was used as the working liquid,the total height of liquid in the tank is 500 mm.Air(ρg=1.225 kg.m-3,μg=1.789×10-5Pa.s)was introduced through a ring sparger (diameter is 210 mm)located below the lower impeller with the gas flow rate Q=0.76 m3.h-1and 1.52 m3.h-1,respectively.There are 20 downward-facing holes(diameter is 2 mm) on the sparger pipe.Fig.4 shows the three impellers used in the experiment.Six-bent-bladed turbine(6BT)as the lower impeller(diameter is Dd=Ts/2),six-inclined-blade down-pumping turbine(6ITD)and six-inclined-blade up-pumping turbine(6ITU)as the upper impeller (diameter is Dd=Du=Ts/2),respectively,i.e.6BT+6ITDand 6BT+6ITU.The back angle (θd) of 6BT blade is 30°,and the inclined angles (θu) of 6ITUand 6ITDblade are all 45°.All of the blade thicknesses are 2 mm.The other dimensions are given in Table 1.

2.2.Experimental results

The distributions of bubbles in stirred tank with two combined dual-impeller under four operating conditions are shown in Fig.5.It can be seen that bubbles distribute mostly in the upper region of stirred tank with large sizes at N=60 r.min-1and Q=0.76 m3.h-1(Fig.5(a)).As Q is increased to 1.52 m3.h-1(Fig.5(b)),the bubbles sizes stirred by 6BT+6ITUimpeller decrease obviously with the broader distribution compared with that by 6BT+6ITDimpeller,but bubbles in bottom region are all less for two combinations.When the speed is increased to 120 r.min-1(Fig.5(c),(d)),it is found that the large bubbles disappear and the even small bubbles appear with the wider distribution.The bubbles tend to be refined with the increase of gas flow rate (Fig.5(d)).

PIV measurement was employed to capture bubbles sizes distributions at different speeds and gas flow rates.The original image of bubbles distributions within axial lengthwise section obtained by PIV was processed into the binary image of black and white with a width of 70 mm and a height of 170 mm,which is shown in Fig.6 at the speed of 60 r.min-1and gas flow rate of 0.76 m3.h-1.The shaft radius in the binary image is shown as 5 mm,while its actual radius is 12.5 mm.The scale may be established based on this ratio of the shaft radius to obtain the real bubbles areas through converting bubbles sizes taken by PIV measurement.

It can be seen from Fig.6 that the bubbles are not circular in shape,so it is necessary to convert the cross-sectional area A of the bubbles into the equivalent circle area.The morphological processing of the captured bubble image was carried out with MATLAB software to obtain the cross-sectional area A of bubbles,as shown in Fig.7.The bubble diameter can be presented according to the equivalent diameter d [30],as defined by Eq.(1).

The diameters of bubbles in five regions from A to E after image processing are shown in Table 2 and Fig.8.It is found from the image post-process results that the bubble distribution ranges from 1 mm to 7 mm.

Table 1Dimensions of the stirred equipment

Table 2Diameters of bubbles in A–E regions for 6BT+6ITU and 6BT+6ITD

Table 3Bubble diameter of each size group used in population balance equations

3.Numerical Simulation

3.1.Computational grid

A pre-processor (Gambit 2.4.6) was used to discretize the flow domain with tetrahedral mesh,and the domain was segmented into the stator zone and rotor zone,as shown in Fig.9.The sufficient amount of mesh nodes was created in the rotor zone including the gas sparger in order to have a much fined mesh.The original 3D mesh had about 345,100 cells.To verify the grid independence[31],the number of cells was increased to about 667,200 and 856,970,respectively.These increases changed the velocity and power number in rotor region by more than 3%.Finally,the number of cells was further increased to 1,104,258,the powernumber and velocity changed less than 3%.Therefore,1,104,258 cells were employed to numerical simulations.

Fig.2.Schematic of PIV experimental apparatus.

3.2.Simulation methods

Fluent 17.0 was applied to simulate the gas–liquid flow in stirred tank.The speed was set between 60 r.min-1and 120 r.min-1in CFD calculation,the corresponding Reynolds numbers range from 138,127 to 276,254.The standard k-ε turbulence model was employed to simulate gas–liquid turbulent flow in stirred tank[27,32,33].The gas–liquid mixing was simulated using the Eulerian-Eulerian multiphase model.

Fig.3.Schematic diagram of gas–liquid stirred tank.

Mass conservation equation for phase i,as given by Eq.(2).The tap water(liquid phase)and air(gas phase)were set as primary phase and secondary phase,respectively.The two phases are assumed to share the space of stirred tank according to their proportion of volume,as given by Eq.(3).

Momentum conservation equation for phase i is computed by Eq.(4)

where Fg,lis interphase force,p′is modified pressure,given by Eq.(5).

For the continuous phase (liquid phase),the effective viscosity is calculated as

Fig.4.Three impellers used in experiment:(a) 6BT,(b) 6ITU,(c) 6ITD.

Fig.5.Bubbles distribution at different conditions.(a) N=60 r.min-1,Q=0.76 m3.h-1 (b) N=60 r.min-1,Q=1.52 m3.h-1 (c) N=120 r.min-1,Q=0.76 m3.h-1 (d) N=120 r.min-1,Q=1.52 m3.h-1.

where μlis the liquid viscosity,μtland μtgrepresent the liquid phase turbulence viscosity and the gas phase induced turbulence viscosity,as given by Eq.(7) and Eq.(8),respectively.

Fig.6.Binary image of bubble distribution.

Fig.7.Post-processing image of bubbles.

where ρ is density,klis turbulent kinetic energy,Cμ,b is empirical parameters with a value of 0.6 [34],and u is velocity vector.

3.3.The drag force

Only the drag contribution is considered as the interphase force,while other forces(lift and virtual mass force)have been neglected because of their less significance in the interaction of phases.The grace drag model was used to solve the interphase force owing to its well calculation accuracy even in the case of different bubble shapes [35,36].The drag force can be obtained by Eq.(9).

where dbis the diameter of the bubble,CDis the drag coefficient.

Different shape flow regimes have relative calculation formulas[37]:

3.4.Simulation procession

The Multiple Reference Frame (MRF) method was firstly adopted to carry out steady-state calculation for only liquid phase flow at the beginning until it converged,and then two-phase flow was conducted as the transient process by sliding mesh (SM)method after air introduction.The time step size of SM was set as 0.001 s and the computational time was 126 hours by the Dell workstation with 2.67 GHz processors and 32 GB RAM.Finally,Population Balance Model(PBM)was coupled to solve the distributions of the bubbles sizes and it has been widely applied to predict the bubble size and its distribution in multiphase systems[34,38,39].The Population Balance Equation (PBE) is as follows:

Fig.8.Distributions of mean bubbles diameters in A–E regions.

Fig.9.Computational grids of gas–liquid stirred equipment.

where n(V,t)is bubbles number density function,a(V,V′)is bubble aggregation rate,g (V′) is bubble breakage frequency,β (V|V′) is probability density function of bubble breakage.Considering the coalescence and breakage phenomena of the bubbles in the gas–liquid stirred tank,the Luo coalescence and Luo &Svendsen breakage model were adopted in calculation owing to its successful application in the gas–liquid two-phase mixing [40–42].The computational scheme is illustrated in Fig.10.The calculations of the flow field and gas holdup were completed first by CFD,and then adding PBM model to CFD.The distributions of bubbles sizes and gas–liquid interfacial areas were calculated by the CFD-PBM coupled model.

Fig.10.Computational scheme of the coupled CFD-PBM model.

Fig.11.Comparison of CFD and experimental power values at different speeds.

3.5.Simulation details

The top surface was modeled as degassing boundary,and the sparger holes were set as the velocity inlet for gas.The momentum equation was discretized by the second-order upwind scheme.The scaled residuals for all variables were less than 1×10-4as convergence condition.The Phase coupled SIMPLE algorithm was performed to couple pressure and velocity terms.The velocities of the liquid and gas phase,as well as gas volume fraction at initial state were all set as zero.

The homogenous MUSIG model was used in simulation.The polydispersed gas phase was divided into 16 sub-size groups ranging from 0.5 to 8 mm based on the experimental results according to the criteria of equal diameter discretization[34,43].The diameter of group i is calculated from equation as where dminand dmaxcorresponds to the minimum and maximum diameter of the gas phase,respectively,and n stands for number of bubble groups.

Fig.12.Distribution of Sauter mean diameters of bubbles from A to E regions.

Fig.13.Flow field structure within axial lengthwise section at N=120 r.min-1 and Q=0.

Fig.14.Flow field structure within axial lengthwise section at N=120 r.min-1 and Q=0.76 m3.h-1.

The diameter distribution of each size group used in the population balance equations was shown in Table 3.The value of 2.15 mm is taken as the bubble size at sparger inlet in numerical simulation,which is in accordance with the inlet diameter of gas disperser in the experiment.

4.Results and Discussion

4.1.Validation of computational model

The power consumption P can be solved out by using Eq.(16).

where N is the speed and M is the torque of the shaft and impeller,which can be obtained by numerical calculation.

The CFD values of power consumption with 6BT+6ITUcombination impeller were compared with experimental data at speeds from 30 to 180 r.min-1and gas flow rate of 0.76 m3.h-1,as shown in Fig.11.It can be seen that CFD simulation results are close to the experimental data with the consistent trend of change and the maximum deviation is only 6.3%.

Fig.12 showed the distribution of Sauter mean diameters of bubbles from A to E regions in stirred tank with two combined impeller at the speed of 120 r.min-1and gas flow rate of 0.76 m3.h-1.The results show that the calculated bubble sizes are always well agreement with the experimental values,and the deviations are reasonable.

Fig.15.Contours of gas holdup within axial lengthwise section:(a) 6BT+6ITU,(b) 6BT+6ITD.

Fig.16.Axial distributions of gas holdup for two combined impellers.

Fig.17.Distribution contours of bubble size within axial lengthwise section:(a) 6BT+6ITU,(b) 6BT+6ITD.

Fig.18.The change curves of the bubbles diameters along the radial direction at different heights.

Fig.19.Distribution proportions of bubbles sizes in stirred tank.

It can be inferred that the CFD-PBM coupled model is validated based on the experiments of the power consumption and distribution of bubbles sizes,which indicates the accuracy and reliability of the CFD-PBM coupled model established.

4.2.Flow field structures

The velocity vector plots within axial lengthwise section at the speed of 120 r.min-1and gas flow rates of 0 and 0.76 m3.h-1are shown in Fig.13 and Fig.14,respectively.It can be seen that the overall circulation of fluid is obviously enhanced at Q=0 for the combined impeller of 6BT+6ITDcompared with 6BT+6ITU,in particular,a strong connecting flow is formed between the dualimpeller,which may greatly enhance the axial flow of the fluid.Whereas the flow field structure changed with the introduction of gas.For 6BT+6ITUcombination impeller,the velocities between the dual-impeller are increased significantly,and the overall circulation of fluid is also strengthened at Q=0.76 m3.h-1,but the change is reversed for 6BT+6ITDat Q=0.76 m3.h-1.So the introduction of gas can contribute to the overall flow of fluid in stirred tank with 6BT+6ITUcombination impeller.

Fig.20.Gas–liquid interfacial area within axial lengthwise section:(a) 6BT+6ITU,(b) 6BT+6ITD.

Fig.21.The change curves of the interfacial area along the radial direction at different heights.

4.3.Gas holdup distribution

The gas holdup is usually used as an important index to characterize the gas–liquid dispersion in a stirred tank.The contours of gas holdup are shown in Fig.15 for two combined impellers at N=120 r.min-1and Q=0.76 m3.h-1.The distributions of gas holdup along axial direction at different gas flow rates and speeds are shown in Fig.16.It can be seen that the overall distributions of gas holdup are similar at the same working conditions.In the vicinity of the upper impeller,the values of gas holdup stirred by 6BT+6ITDcombination impeller are generally larger than that by 6BT+6ITUwith a large amount of gas accumulated near the upper impeller and shaft,while the distributions of gas holdup stirred by 6BT+6ITUare more uniform in this region.In the vicinity of the lower impeller,the gas holdup produced by 6BT+6ITUhas a broader distribution.This is mainly because that the upward flow pattern produced by 6ITUimpeller is in the same direction as the gas rises,which can strengthen the gas dispersion.Moreover,the overall gas holdup are more obviously increased with the speed from 60 to 120 r.min-1compared with the increasing gas flow rate.

4.4.Bubble size and distribution

Bubble size plays an important role for the evaluation of gas–liquid interfacial area and momentum exchange between phases and is a key index to assess the mixing performance of the gas–liquid stirred tank.Fig.17 shows the distribution contours of bubbles sizes at N=120 r.min-1and Q=0.76 m3.h-1.The change curves of the bubbles diameters along the radial direction at different heights are showed in Fig.18.It can be seen that the larger sizes of bubbles produced by 6BT+6ITDcombination impeller are mainly concentrated near the upper impeller and upper area of shaft,while the bubbles sizes produced by 6BT+6ITUare smaller and have a more uniform distribution in the corresponding region.Besides,in the area below the upper impeller,it is also found that the large bubbles produced by 6BT+6ITUaccount for a low proportion with the more uniform distribution of bubbles sizes compared with that by 6BT+6ITD.This is mainly because that the flow field structure produced by 6BT+6ITUis conducive to the dispersion of bubbles.

Fig.19 shows the distribution proportions of bubbles sizes in stirred tank with two combined impellers.It can be observed that the 4 mm bubbles account for the largest proportion of 33% at Q=0.76 m3.h-1and 22% at Q=1.52 m3.h-1for 6BT+6ITU,while the bubbles of 4.7 mm and 5.5 mm are the largest proportion for 6BT+6ITD,i.e.the corresponding proportions are 25% at Q=0.76 m3.h-1and 22%at Q=1.52 m3.h-1,respectively.At the meantime,the proportions of other bubbles sizes (0.6–3.6 mm and 5.7–7.4 mm)are close for two combinations.So 6BT+6ITUcombination impeller can suppress the large bubbles and reduce bubbles sizes,effectively,as well as promote gas dispersion.

4.5.Gas–liquid interfacial areas

The gas–liquid interfacial area can be used as an important index to characterize the mass transfer rate.It can be expressed by the following formula

where αgis the gas holdup,d32is Sauter mean diameter.

The contours of gas–liquid interfacial area at N=120 r.min-1and Q=0.76 m3.h-1are showed in Fig.20.The change curves of the interfacial area along the radial direction at different heights are showed in Fig.21.It can be observed that there is gas accumulation near upper area of shaft and the upper impeller region for 6BT+6ITDcombination impeller,which is not conducive to gas–liquid mass transfer,while the gas–liquid interfacial area produced by 6BT+6ITUhas a broader and uniform distribution without gas accumulation in the whole tank.

5.Conclusions

The coupling algorithm of CFD and PBM model was carried out to investigate gas–liquid two-phase flow characteristics in stirred tank with two combined dual-impellers,i.e.6BT+6ITUand 6BT+6ITD.The effects of speed and gas flow rate on the flow field structure,gas hold-up,distribution of bubbles sizes and gas–liquid interfacial area were analyzed comparatively.The conclusions are obtained as following:

(1) The introduction of gas will change the flow field structure and can facilitate the mixing of fluid in stirred tank with 6BT+6ITUcombination impeller.

(2) 6BT+6ITUcombination impeller can make more uniform bubbles sizes with a better bubble dispersion performance compared with 6BT+6ITDcombination impeller.Moreover,the distributions of gas holdup and the gas–liquid interfacial area produced by 6BT+6ITUare more even and wider.

(3) The gas holdup will be improved obviously with the increasing speed compared with the gas flow rate,which indicates that the high speed is more beneficial to the increase of gas holdup.

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

This work was supported by the National Natural Science Foundation of China (52176040) and Shandong Provincial Natural Science Foundation of China (ZR2018LE015).

Nomenclature

A bubble cross-sectional area,m2

a gas–liquid interfacial area,m2

b clearance from the lower impeller to the bottom of tank,m

bdimpeller width,m

c clearance from the sparger to the bottom of tank,m

D the spacing between the two impellers,m

Ddlower impeller shaft diameter,m

Duupper impeller shaft diameter,m

d bubble equivalent diameter,m

dzstirring shaft diameter,m

d32Saute mean diameter,m

Eo E?tv?s number

FDinterphase force,Pa

g gravitational force,m.s-2

hbbaffle height,m

klliquid turbulence energy,m2.s-2

N rotational speed of impeller,r.min-1

P power consumption,W

p′modified pressure,Pa

Q gas flow rate,m3.h-1

Re Reynolds number

Tsvessel diameter,m

u velocity vector,m.s-1

w baffle width,m

α phase volume fraction,gas holdup

μ dynamic viscosity,Pa.s

μeff,leffective viscosity of liquid phase,Pa.s

ρ density,kg.m-3

Subscripts

cap-bubble cap-shape bubble

ellipse ellipse-shape bubble

g gas

l liquid

sphere spherical shape bubble

主站蜘蛛池模板: 91精品国产自产91精品资源| 亚洲天堂网在线观看视频| 91人妻在线视频| 国产一区二区三区精品欧美日韩| 久久精品这里只有国产中文精品| 精品伊人久久久香线蕉| 在线看片免费人成视久网下载| 日韩欧美国产成人| 亚洲a免费| 免费播放毛片| 日韩美毛片| 亚洲成人在线网| 在线播放国产一区| 蜜臀av性久久久久蜜臀aⅴ麻豆| 亚洲A∨无码精品午夜在线观看| 国产精品免费电影| 久视频免费精品6| 国产农村精品一级毛片视频| 1级黄色毛片| 午夜无码一区二区三区在线app| 青青草原国产av福利网站| 亚洲精品成人片在线观看| 免费观看欧美性一级| 97人人模人人爽人人喊小说| 免费无码又爽又刺激高| 国产农村1级毛片| 亚洲无码精品在线播放| 亚洲福利一区二区三区| 日本国产在线| 9久久伊人精品综合| 青草视频网站在线观看| 亚洲欧美综合在线观看| 国产呦精品一区二区三区下载 | 尤物在线观看乱码| 午夜无码一区二区三区| 91成人免费观看在线观看| 99re视频在线| 国产青青操| 国产情侣一区二区三区| 2020亚洲精品无码| 久久久久久午夜精品| 亚洲精品无码不卡在线播放| 亚洲一区免费看| 国内精品视频在线| 精品剧情v国产在线观看| 免费又爽又刺激高潮网址 | 香蕉久人久人青草青草| 91尤物国产尤物福利在线| 狠狠ⅴ日韩v欧美v天堂| 亚洲天堂网视频| 91视频99| 午夜视频免费一区二区在线看| 亚洲AⅤ永久无码精品毛片| 国产18在线| 综合网久久| 天堂亚洲网| 欧美日本在线| 日本道综合一本久久久88| 天天操精品| 伊人久热这里只有精品视频99| 天堂av综合网| 久久亚洲国产一区二区| 久久这里只有精品免费| 91精品国产91久久久久久三级| 国产精品美女网站| 欧美日韩另类在线| 久久国产黑丝袜视频| 日韩在线1| 99资源在线| 精品国产成人高清在线| 亚洲成人动漫在线观看| 4虎影视国产在线观看精品| 女人18毛片久久| 久久综合九色综合97婷婷| 亚洲AⅤ无码日韩AV无码网站| 亚洲乱码在线视频| 伊人久久婷婷五月综合97色| 国产白浆一区二区三区视频在线 | 97人妻精品专区久久久久| 园内精品自拍视频在线播放| 日韩午夜福利在线观看| 国产玖玖玖精品视频|