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

Numerical simulation of flow and heat transfer of n-decane in sub-millimeter spiral tube at supercritical pressure

2023-11-12 02:01:04JiahaoXingHuaizhiHanRuitianYuWenLuo

Jiahao Xing,Huaizhi Han*,Ruitian Yu,Wen Luo

School of Chemical Engineering,Sichuan University,Chengdu 610065,China

Engineering Research Center of Combustion and Cooling for Aerospace Power,Ministry of Education,Sichuan University,Chengdu 610065,China

Keywords:Sub-millimeter spiral tube Supercritical pressure Numerical simulation Heat transfer performance

ABSTRACT The flow and heat transfer characteristics of n-decane in the sub-millimeter spiral tube(SMST)at supercritical pressure (p=3 MPa) are studied by the RNG k -ε numerical model in this paper.The effects of various Reynolds numbers(Re)and structural parameters pitch(s)and spiral diameter(D)are analyzed.Results indicate that the average Nusselt numberand friction factor(ˉf)increase with an increase in Re,and decrease with an increase in D/d (tube diameter).In terms of the structural parameter s/d,it is found that as s/d increases,the Nu and ˉf first increase,and then decrease.and the critical structural parameter is s/d=4.Compared with the straight tube,the SMST can improveby 34.8%at best,while it can improve ˉf by 102.1% at most.In addition,a comprehensive heat transfer coefficient is applied to analyze the thermodynamic properties of SMST.With the optimal structural parameters of D/d=6 and s/d=4,the comprehensive heat transfer factor of supercritical pressure hydrocarbon fuel in the SMST can reach 1.074.At last,correlations of the average Nusselt number and friction factor are developed to predict the flow and heat transfer of n-decane at supercritical pressure.

1.Introduction

As modern science and technology have advanced,the pressure ratio of the aero-engine and the turbine inlet temperature are continually increasing and the cooling system of turbine engines is greatly challenged.In order to meet this stringent cooling demand at present,an internal cooling technology named CCA(cooled cooling air) has been proposed.In this technology,the supercritical pressure hydrocarbon fuel is used to cool the compressed air [1].Due to the high density and heat sink properties of hydrocarbon fuels [2],the technology can greatly reduce the air temperature,significantly improving the air cooling efficiency.Besides,it can increase the fuel temperature at the inlet of the combustion chamber,which promotes atomization and combustion compared with using conventional coolant [3,4].In addition,continued research into new cooling structures is necessary.And various structures have been studied by scholars [5-7].Due to their compact structure,high efficiency,and ease of integration,spiral tube heat exchangers are widely used as air coolers.At the same time,in order to further increase the compactness of the heat exchanger,the size of the heat exchange tube is gradually reduced.Besides,for supercritical pressure fluids,tube size is one of the important parameters affecting its flow heat transfer characteristics due to the existence of flow acceleration effect and buoyancy effect [8].As a consequence,to efficiently design the sub-millimeter spiral tube(SMST)with supercritical pressure fluid,it is necessary to recognize their special flow and heat transfer performance.

Under supercritical pressure,the physical properties of hydrocarbon fuels will change greatly with the increase of temperature,making its heat transfer characteristics very complex and unique compared with the fluid with constant physical properties.Fig.1 shows the changes in the physical properties ofn-decane under supercritical conditions,which is the main component of hydrocarbon fuel.There have been many studies on the flow and heat transfer performance of supercritical fluids.Lietal.[9] carried out a numerical study on the convective heat transfer characteristics of Chinese aviation kerosene RP-3 under supercritical pressure by using the finite volume method.Xuetal.[10]found that the ribbed tube has a good promotion effect on the heat transfer of hydrocarbon fuel under supercritical pressure,and the tube wall temperature can be reduced by about 300 K in the inlet region.Fengetal.[11] inserted a twisted band in a circular tube and studied the convective heat transfer of hydrocarbon fuels in the tube under supercritical pressure.Yanetal.[12]conducted experiments on the flow and heat transfer ofn-decane in a vertical tube under supercritical pressure.By slowly increasing the heat flow,different phases of flow are observed.Jiangetal.[13] found that supercritical pressure hydrocarbon fuels with thermal cracking benefit from an additional heat sink near the wall,resulting in enhanced heat transfer.But at the same time,the convective heat transfer is weakened due to the increase in thermal resistance due to coking and bubble layers.Wenetal.[14] conducted an experimental study on the convective heat transfer characteristics of aviation kerosene RP-3 in a vertical spiral tube under supercritical pressure.And two kinds of Nusselt correlations of aviation kerosene RP-3 in spiral tubes under supercritical pressure are established.

Fig.1.Physical properties of n-decane.

The spiral tube structure was widely used in the 1890s.And spiral heat exchangers play an important role in steam generators,chemical plants,refrigeration systems,ships,and nuclear reactors[15,16].As shown in Fig.2,the main application of this device in turbine engines is in the external duct,which is the blue part in the figure.As early as the beginning of the 20th century,Dean researched fluid flow in spiral tubes and predicted the formation of secondary flows [17,18].Besides,a numerical study of incompressible flow and heat transfer in a three-dimensional fixed curved pipe by Nobari and Amani [19].Ciofaloetal.[20] studied the effects of gravity and centrifugal buoyancy on laminar flow and heat transfer in the spiral tube by numerical simulation.The experimental results of Prabhanjanetal.[21]on spiral and straight tubes of equivalent length show that the heat transfer coefficient of the spiral tube is larger.Naphonetal.[22] investigated the curvature ratio of the horizontal spiral coil through experiments and numerical simulations.The results show that both heat transfer and pressure drop in the spiral tube are greater than in the straight tube.Keetal.[23]used numerical simulation to analyze the influence of structural parameters such as cone angle,section shape,and pitch on the heat transfer characteristics of the conical spiral tube.Sasmitoetal.[24] studied the heat transfer performance of in-plane spiral tubes with different cross-sectional shapes under constant wall temperature and wall heat flux density.Moreover,Kurnia and Sasmito [25] found that the heat transfer performance of square spiral tubes with the same straightened length increases gradually with the decrease of the radius of curvature.Zhaietal.[26] conducted experiments on the flow and heat transfer of nanofluids in spiral tubes,and the influence of the pitch and rotation angle is analyzed.Lietal.[27] established a numerical model of a spiral microtube reactor and studied the effects of pipe diameter,pitch,and spiral diameter on its heat transfer performance.

Fig.2.CCA heat exchanger in turbine engines.

Meanwhile,attention has been paid to the study of submillimeter tubes.It has been extensively studied by scholars[28-31].Wenetal.[32]conducted experiments on the heat transfer of supercritical aviation hydrocarbon fuel RP-3 in horizontal miniature tubes and found that even small round tubes have a greater impact on the buoyancy effect.Zhouetal.[33] found that the flow instability of supercritical hydrocarbon fuel during small-scale tube cooling can be improved by reducing the mass flow rate and using a fuel whose density changes more slowly with temperature.Fuetal.[34] studied the convective heat transfer of supercritical pressure hydrocarbon fuel RP-3 in vertical submillimeter tubes with inner diameters ranging from 0.5 to 2 mm.While Puetal.[35] studied the convective heat transfer and flow resistance characteristics of supercritical pressure hydrocarbon fuels in horizontal rectangular mini-tubes.Huaetal.[36] found through numerical simulation that when the inlet velocity is greater than 10 m.s-1,the supercritical heat transfer coefficient ofn-heptane in horizontal micropipes can be predicted by the supercritical heat transfer expression of CO2.Leietal.[37] carried out numerical simulations on the cooling heat transfer and pressure drop of supercritical CO2in wavy microchannels.Zhaoetal.[8]found that the heat transfer of supercritical CO2in vertical pipes with inner diameters of 0.27 mm is degraded due to the flow acceleration effect,while Yangetal.[38] showed that the heat transfer and pressure drop performance of supercritical CO2in a microchannel of about 0.5 mm is not much different from that in the normal case.Khalesi and Sarunac [39] also obtained similar conclusions from the numerical analysis of supercritical CO2flow and conjugated heat transfer in square microchannels.Jajjaetal.[40] considered the need for systematical analysis of microscale heat transfer under non-uniform heating conditions.

References listed above focus primarily on the flow and heat transfer of supercritical fluids in conventional-scale spiral tubes and sub-millimeter straight tubes.However,there are few studies on the flow and heat transfer performance of supercritical pressure fluids in the sub-millimeter spiral tubes (SMSTs).Therefore,the purpose of this paper is to identify the effects of structural parameters,such as Reynolds number (Re),the ratio of pitch to tube diameter (s/d),and the ratio of spiral diameter to tube diameter(D/d) on supercritical fluid flow and heat transfer in the SMSTs.It can provide a significant design consideration for the SMST,supporting the miniaturization and compactness of aerospace heat transfer technology.

2.Physical Model and Numerical Method

2.1.Physical model and grid generation of the SMST

As shown in Fig.3,SMSTs with different structural parameters are constructed in this study.The influence ofs/d,D/d,andReon the flow and heat transfer characteristics of the tube are studied.The SMST of different structures have the same length(L=120 mm),and their inner diameters (d) are fixed at 0.5 mm.Table 1 summarizes all influence parameters.

Table 1 Influence parameters of the SMST

Fig.3.Physical model of the SMST (unit: mm).

Modeling and meshing are then performed using the ICEM software.A high-quality grid is obtained by controlling the aspect ratio of the grid within 200 and dividing the mainstream area grid using the‘‘O-Block”method.Lastly,the grid quality is above 0.5,which is high quality and good for computation.Additionally,due to the large velocity gradients and temperature gradients,encryption processing is performed in the near-wall region as illustrated in Fig.4.And the wall enhancement function is used to calculate the flow and heat transfer in the boundary layer,which requires they+at the wall to be close to 1.

Fig.4.The grid division of the SMST.

2.2.Governing equations

The governing equations are given as follows:

(1) Continuity equation:

(2) Momentum equation:

(3) Energy equation:

In this paper,thek-ε RNG turbulence model is selected as the numerical simulation method.For thek-ε RNG turbulence model,kand ε are two basic unknown quantities,and the corresponding transport equations are whereGkis the turbulent kinetic energy generated by the average velocity gradient.αk,αεare the effective Prandtl numbers of thekequation and ε equation,respectively.AndC1ε,C2ε are constant.

2.3.Boundary condition

ANSYS FLUENT 19.0 is used for numerical simulation.The pressure-based and steady-state solver solves the governing equations,and the COUPLED algorithm is chosen for pressure-velocity coupling.The discrete schemes of pressure,momentum,energy,turbulent kinetic energy,and specific dissipation rate are defined as second-order upwind.Except that the residual error of the continuity equation is set to 10-3,the other residual error of the equation is set to 10-6.

The following are detailed boundary conditions:

(1) Inlet Reynolds numbers range from 3000 to 5000 in increments of 500,with the fluid fully developed.

(2) The inlet temperature isTin=300 K.

(3) Pressure outlet is imposed at the outletPout=3 MPa.

(4) Constant heat flux is imposed at the channel wallqw=0.8 MW.m-2.

(5) Non-slipping boundary condition at the channel walluw=0.

2.4.Grid independence and numerical model validation

In this paper,the average convective heat transfer coefficient of the tubeis calculated by:

The expressions of the average Nusselt numberand the friction factorare:

The comprehensive heat transfer coefficient η is adopted to evaluate the comprehensive heat transfer performance for the SMST.When η is greater than one,it indicates that the comprehensive heat transfer performance of the SMST is better than the straight tube.The expression of η is:

In order to improve the computing efficiency while guaranteeing the computing accuracy,grid independence of the current simulations has been conducted.Considering the large fluid velocity and temperature gradient in the boundary layer region,the mesh thickness of the first layer is changed to control the near wall meshing number.The number of grids increases gradually as the thickness of the first layer decreases,resulting in a gradually decrease ofIt can be seen from Fig.5,when the thickness of the first layer reduces to 0.0005 mm,the results are basically stable,and there are about 1.45 million grids at this time.Therefore,the first layer of mesh thickness is chosen as 0.0005 mm for the simulation.Table 2 provides a detailed description of the grid parameters.The value ofy+is determined by the thickness of the first layer of the grid.By verifying the grid independence,we determined the thickness of the first grid layer.At this thickness,y+at the wall is approximately 1.All of the final grids chosen have ay+below 1.

Fig.5.Grid independence.

The numerical simulation results are compared with the experimental results of Wenetal.[14] to verify the accuracy of the numerical simulation method used in this paper.In the experiments and simulations,the spiral tube has a total length of 1500 mm in length,with a spiral diameter of 20 mm,a pitch of 10 mm,and a constant heat flux of 240 kW at the tube wall.Fig.6 compares the average Nusselt number calculated by different numerical models with the experimental results.It shows that thek-ε RNG turbulence model is the most consistent model with the experimental results,with an error of within 12%.Hence,thek-ε RNG turbulence model used in the paper is accurate and feasible.

Fig.6.Model validation.

Fig.7.Effects of D/d on the heat transfer: (a) (b)

Fig.8.Effects of s/d on the heat transfer: (a) (b)

3.Results and Discussion

3.1.Heat transfer performance for SMST

3.1.1.EffectsofD/donheattransfer

3.1.2.Effectsofs/Donheattransfer

3.2.Flow performance for SMST

3.2.1.EffectsofD/dontheflowperformance

Fig.9 shows the effects ofD/don theof the SMST with variousRe.It is found from Fig.9(a) that the resistance performance of the SMST gradually gets better asD/dincreases,and the trend of getting better is gradually slowing down.Besides,it presents that the maximumatD/d=4 is 45%higher than the minimum atD/d=16.From Fig.9(b),it can be seen that with an increase inalso declines gradually.Meanwhile,the SMST has a maximumof 2.02 times of a straight tube whenD/d=4.

3.2.2.Effectsofs/dontheresistanceperformance

Effects ofs/don the flow performance are researched,too.Fig.10 presents the effects ofs/don theof SMST with differentRe.The results show that the resistance performance deteriorates in the first small range with increasings/d,and then gradually improves ass/dincreases.The maximum and minimum values ofin Fig.10(a)differ by at most 12.8%.The laws ofandare similar.Besides,thefor SMST can reach 1.72 times at most that of the straight tube whens/d=4 as shown in Fig.10(b).

Fig.10.Effects of s/d on the resistance performance: (a) (b)

3.3.Comprehensive heat transfer performance evaluation

Based on the above analysis,it is concluded that the structure of SMST not only results in an enhanced heat transfer performance,but also results in a higher resistance loss.In order to comprehensively evaluate the thermo-hydraulic performance of the SMST,a comprehensive heat transfer factor(η)based on Eq.(10)is applied in this study.Fig.11 gives the effects ofs/dandD/don the comprehensive heat transfer factor respectively.The results indicate that η decreases with the increasingRe.As well as this,the SMST is able to achieve the best comprehensive heat transfer performance by optimizing the structural parameterss/d(s/d=4) andD/d(D/d=6).It reveals spiral structures can effectively strengthen heat transfer before the critical structural parameters (s/d<4,D/d<6)without causing excessive resistance losses.However,after the critical structural parameters (s/d>4,D/d>6),it must consume much more resistance coefficient to achieve the heat transfer enhancement.According to the results,the maximum comprehensive heat transfer coefficient of the SMST can reach 1.074 at best whens/d=4 andD/d=6.

Fig.11.Effects of structural parameters on the flow comprehensive heat transfer performance: (a) s/d,(b) D/d.

3.4.Flow and heat transfer mechanism for SMST

3.4.1.AnalysisofD/d

Fig.12 contains the velocity contours of the outlet,as well as the horizontal(a)and vertical(b)velocity distribution of the center section of the SMST with differentD/dwhens/d=8.

Fig.12.Outlet velocity contours and the centerline velocity distributions with different D/d: (a) velocity contours,(b) horizontal velocity distribution,(c) vertical velocity distribution.

Fig.12(a)illustrates that the fluid in the SMST with smallerD/dhas a higher velocity close to the outside of the SMST.Furthermore,it can be speculated that the velocity of the fluid near the inside of the SMST with largerD/dis greater due to the conservation of momentum,which can be proved by Fig.12(b).And it leads to a larger velocity gradient close to the outside and a smaller velocity gradient near the inside of the SMST with greaterD/d.

In addition,it can be observed from Fig.12(c) that the fluid velocity in the upper and lower parts of the SMST with smallerD/dis larger and closer to the wall than in the SMST with largerD/d.According to Fig.12(c),this is reflected in a more significant M-shaped velocity distribution,which indicates a greater velocity gradient on the wall.

The velocity vector diagrams of differentD/dare shown in Fig.13.It is evident that fluid in the SMST with smallerD/dhas a greater flow velocity towards the outer side wall,and the velocity gradient is larger at the wall.The results indicate that the SMST with smallerD/dexhibits a greater secondary flow effect.

Fig.13.Outlet velocity vector diagrams with different D/d.

In summary,it can be concluded that centrifugal force has a greater effect on the SMST with smallerD/d.In addition to increasing the velocity gradient at the wall,it also enhances the secondary flow effect and promotes the mixing of fluids.Both of the above will degrade the resistance performance.

Fig.14 shows the temperature contours of the fluid at the outlet of the SMST with differentD/d,as well as the horizontal and vertical temperature distributions of the center section whens/d=8.

Fig.14.Outlet temperature contours and centerline temperature distributions with different D/d: (a) temperature contours,(b) horizontal temperature distribution,(c)vertical temperature distribution.

In Fig.14(a),it is evident that the temperature of the fluid near the outside of the SMST increases gradually with increasingD/d.And it can be inferred from the law of energy conservation that the temperature of the fluid near the inner side of the SMST gradually decreases asD/dincreases,as shown in Fig.14(b).Further,Fig.14(b) indicates that the fluid in the SMST with smallerD/dhas a thicker thermal boundary layer near the inside of the SMST.Conversely,the thermal boundary layer is thinner near the outside of the SMST whenD/dis smaller.

Moreover,a comprehensive analysis of Fig.14(a) and Fig.14(c)indicates that the fluid in the SMST with smallerD/dhas a lower temperature and thinner thermal boundary layer thickness in the upper and lower parts of the SMST,which appears as a more pronounced W-shaped velocity distribution in Fig.14(c).

To sum up the above,the SMST with a smallerD/dhas a thicker thermal boundary close to the inside of the SMST,while it is thinner elsewhere.Consequently,the temperature of the fluid is reduced in most areas of the SMST with smallerD/d,which results in better heat transfer performance.

In Fig.15,the turbulent kinetic energy contours of the fluid at the outlet of SMST are shown with differentD/d.AsD/ddecreases,the turbulent kinetic energy of the fluid near the inside of the SMST gradually decreases,indicating that the fluid is more intensely mixed,thereby improving heat transfer.

Fig.15.Outlet turbulent kinetic energy contours with different D/d.

3.4.2.Analysisofs/d

Fig.16 contains the velocity contours at the outlet,as well as the horizontal and vertical velocity distributions of the center section of the SMST with differents/dwhenD/d=8.

The difference in Fig.16(a)is difficult to discern with the naked eye.Also,the difference between the curves in Fig.16(b) is minimal.Nevertheless,it can be observed that the fluid velocity in the SMST withs/d=4 is slightly higher than that in the SMST withs/d=2 ands/d=16 near the upper wall of Fig.16(c).Due to this,the velocity gradient at the SMST wall withs/d=4 is slightly larger,resulting in greater resistance.

The velocity vector diagrams of differents/dare shown in Fig.17.Whens/dis larger (s/d=14),the velocity vector in the mainstream region does not point outward,but is bent instead obviously.Therefore,it reveals that the effect of centrifugal force is insufficient to fully develop the secondary flow in the SMST with larges/d.Furthermore,whens/dis larger (s/d=14),the velocity flowing to the outer wall of the pipe is smaller,indicating a weaker secondary flow.

Fig.17.Outlet velocity vector diagram with different s/d.

Overall,the parameters/dhas little effect on the fluid flow performance of the SMST.In the SMST,flow resistance increases ass/ddecreases.However,whens/d<4,frictional resistance decreases because of a reduction in velocity in the vertical direction.This phenomenon is known as the velocity effect [41].

Fig.18 shows the temperature contours of the fluid at the outlet of the SMST with differents/d,as well as the horizontal and vertical temperature distributions of the center section whenD/d=8.

Fig.18.Outlet temperature and centerline temperature distributions with different s/d: (a) temperature contours,(b) horizontal temperature distribution,(c) vertical temperature distribution.

In accordance with Fig.18(a),there is almost no difference between the temperature contours.As shown in the partially enlarged view of Fig.18(b),the SMST withs/d=4 has a higher temperature and a thicker thermal boundary layer close to the inside of the SMST than the SMST withs/d=2 ands/d=14.According to the law of conservation of energy,it is possible to speculate that temperatures are lower and the boundary layer is thinner near the exterior of the SMST.In addition,the enlarged view of Fig.18(c)indicates that the SMST withs/d=4 has a lower temperature and a thinner thermal boundary layer in the upper part of the SMST.

In general,we can conclude that by decreasings/d,the thickness of the thermal boundary layer of the fluid in the SMST will be slightly reduced,thereby improving heat transfer efficiency.However,due to the velocity effect,the heat transfer performance will suddenly decrease whens/d<4.In conclusion,the structural parameters/dhas little effect on the heat transfer performance of the SMST.

In Fig.19,the turbulent kinetic energy contours of the fluid at the outlet of the SMST are shown with differentD/d.According to Fig.19,the turbulent kinetic energy of the fluid near the inner side of the tube increases ass/ddecreases.This is because the fluid on the inner side of the pipe mixes and exchanges more vigorously ass/ddecreases,which is beneficial to heat transfer.

Fig.19.Outlet turbulent kinetic energy contours with different s/d.

4.Heat Transfer and Flow Correlations

Based on the numerical simulation data,the correlation ofNuandof the supercritical pressuren-decane in SMST is proposed.Since the temperature difference between the inlet and outlet is within 80 K,the influence ofPris ignored.Thence,the following correlations of the predicted datacan be obtained as follows:

A comparison is made between the data from the predicted model and those from the simulated model.As shown in Fig.20,the error between theis within 4%,while the error between theandis within 7%.It shows the validity of Eqs.(11) and (12).

5.Conclusions

This paper numerically studies the flow and heat transfer characteristics of supercritical hydrocarbon fuels in the SMST with different structural parameters.Some primary conclusions are obtained as follows:

(2) With the optimal structural parameters ofD/d=6 ands/d=4,the comprehensive heat transfer factor of supercritical pressure hydrocarbon fuel in the SMST can reach 1.074.

(3) A secondary flow is formed due to the centrifugal force of the fluid in the SMST,which significantly break the thermal boundary layer near the outside of the tube wall.The heat transfer and flow resistance performance with various structural parameters can be deeply reflected by the boundary layer and vortex distribution.

(4) The correlations of the average Nusselt number and friction factor are developed,considering structural parameters(D/dands/d)andRe,in order to predict the flow and heat transfer ofn-decane under supercritical pressure.

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 gratefully acknowledge the support by the Scientific Research Start-up Funds for introducing Talent in the Sichuan University(20822041C4014)and National Science and Technology Major Project of China (2017-I-0004-0004).

Nomenclature

cpspecific heat capacity,J·kg-1.K-1

Dspiral diameter,m

dinner diameter,m

Eenergy,J.m-3

ffriction factor

Gkturbulent kinetic energy,m2.s-2

hheat transfer coefficient,W.m-2.K-1

Iunit tensor

kthermal conductivity,W.m-1.K-1

Loverall length,m

NuNusselt number

ppressure,MPa

qheat flux,MW.m-2

spitch,m

Ttemperature,K

ReReynolds number

α effective Prandtl numbers

ε uncertainty

η comprehensive heat transfer factor

μ dynamic viscosity,Pa.s

v kinetic viscosity,m2.s

ρ density,kg.m-3

τ stress,Pa

Subscripts

b bulk

in inlet

out outlet

s straight

w wall

主站蜘蛛池模板: 国产精品永久在线| 97在线免费视频| 在线不卡免费视频| 性欧美在线| 欧美啪啪精品| 国产精品极品美女自在线网站| 99久久人妻精品免费二区| 国产精品刺激对白在线| 在线免费看黄的网站| 五月天婷婷网亚洲综合在线| 欧美日韩中文国产| 久久婷婷五月综合色一区二区| 色综合天天综合中文网| 国产第一页屁屁影院| 激情网址在线观看| 国产无吗一区二区三区在线欢| 91偷拍一区| 亚卅精品无码久久毛片乌克兰 | 男人的天堂久久精品激情| 国产二级毛片| 日韩美女福利视频| 国产在线观看精品| 九九九精品成人免费视频7| 国产一区免费在线观看| 伊人激情综合网| 日韩国产黄色网站| 亚洲欧州色色免费AV| 爽爽影院十八禁在线观看| 亚洲成人一区二区三区| www.youjizz.com久久| 国产午夜人做人免费视频中文 | 国产免费高清无需播放器| 久夜色精品国产噜噜| 国产又色又爽又黄| 成人精品在线观看| 久久婷婷色综合老司机| 欧美精品高清| 666精品国产精品亚洲| 午夜免费小视频| AV天堂资源福利在线观看| 国产综合另类小说色区色噜噜| 国产成人精品在线1区| 久操中文在线| 久久久久人妻精品一区三寸蜜桃| 全免费a级毛片免费看不卡| 国产区人妖精品人妖精品视频| 麻豆国产精品一二三在线观看| 欧美国产三级| 在线观看国产精品第一区免费| 亚洲激情区| 囯产av无码片毛片一级| 精品国产成人a在线观看| 国产欧美在线观看精品一区污| 亚洲国产日韩欧美在线| 婷婷五月在线| 2021最新国产精品网站| a在线亚洲男人的天堂试看| 国产欧美日韩综合一区在线播放| 中国一级毛片免费观看| 国产不卡国语在线| 波多野结衣一区二区三区四区视频| 99re免费视频| 国产色网站| a级毛片在线免费观看| 国产在线视频二区| 四虎国产在线观看| 亚洲精品不卡午夜精品| 国产网友愉拍精品| 欧洲成人在线观看| 69av在线| www.av男人.com| 91欧美亚洲国产五月天| 久久久久久久久久国产精品| 国产亚洲视频中文字幕视频 | 91精品久久久无码中文字幕vr| 国产精品久线在线观看| 视频二区亚洲精品| 国产美女视频黄a视频全免费网站| 久热re国产手机在线观看| 国产h视频在线观看视频| AV网站中文| 在线精品视频成人网|