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

    Afterbody aerodynamic optimization design of transport airplane considering wing wake flow

    2015-11-08 06:18:18BaiJunqiangSunZhiweiDongJianhongHuangJiangtao
    關(guān)鍵詞:張彬外形框架

    Bai Junqiang,Sun Zhiwei,Dong Jianhong,Huang Jiangtao

    (1.Northwestern Polytechnical University,Xi’an 710072,China;2.AVIC The First Aircraft Institue,Xi’an 710089,China)

    0 Introduction

    The technique of transport’s afterbody shape design for drag reduction has gained extensive attentionfrom most aircraft manufacturers[1-3].The drag amount of the afterbody is about 1/3 of the total drag of the transport[4].Flow separation will easily occur[1-4]for that the flow field is so complex under the influence of downwash,which will further cut down the safety and economy due to the structural oscillation of the afterbody[5].Hence the research on the design method of transport’s afterbody shape with consideration of flow separation and vortex-induction drag has been paid a lot of attention in many countries,which keep ahead in advanced aeronautical techniques.It is of great value in both theory and engineering to study the influence of afterbody shape on flow characteristics[6],which is the key technique of transport’s body design.Statistical results show that the drag could be reduced by about 0.5%-3%with a good design of the shape,hence remarkable economic benefit will be achieved.

    The upswept angle,fineness,flatness and contraction ratio are several key design parameters in the shape design of the afterbody.Upswept angle mainly determines the flow field separation characteristics which directly relate to the pressure drag and structural oscillation.Most of these research works on optimizing the parameters in fuselage model or 2D configuration without the influence from other parts of the aircraft[6-7].However,when wing is a primary part to affect the flow field with downwash,they will have a great effect on the flow field and the optimized result.Therefore,in order to get the optimized afterbody shape in a whole aircraft model,a more refined optimization frame has to be carried out considering the wing wake flow.

    The main difficulty in afterbody optimization lies on the parametric techniques and optimized frame to reduce the computational cost.In present analysis,a transport’s afterbody shape is optimized by a comprehensive optimization framework that has been established by integration of advanced surrogate model and geometric parameterization.Firstly,F(xiàn)ree Form Deformation(FFD)technique is used for afterbody shape parametric model.Secondly,Kriging surrogate model is for aerodynamic characteristics approximation.Thirdly,infinite interpolation method is for spatial grid deformation.Finally,the afterbody optimization design result shows that the drag coefficient of wingbody configuration is reduced dramatically.

    1 FFD parametric method

    The FFD method,which integrates the deformed modeling method into traditional CAD/CAM system,was proposed in 1986 by Sederberg and Parry of Brigham Young University.It can provide the designers with more freedom in modifying the shape with deformation,and has a good property to ensure smoothness and continuity of the geometry.The algorithm assume that the object has good flexibility,which can be deformed easily under the action of external forces.

    The procedures of FFD method to manipulate theobject’s geometry are as following:

    (1)The first step is creating a deformation tool to define a control frame with control points lying on each of its edges and the local i-j-k as their identity,as shown in Figure 1.

    Fig.1 Creating a deformation frame圖1 構(gòu)建變形框架

    (2)The second step is calculating the original Cartesian coordinates of the research object’s control points and mapping the research object from the Cartesian coordinate system to parameter space.The parameter coordinates of the research object’s geometrical points in the parameter space should be calculated according to the local Cartesian coordinates.

    (3)Modifying the deformation tool.The control points’original coordinates in the tensor product control body can be changed to get the new global coordinates,as shown in Figure 2.

    The new Cartesian coordinates of the research object can be calculated by control points’new global coordinates and the local parameter coordinates of the research object in the constant parameter space,as Figure 3 shows.

    Fig.2 Modifying the deformation frame圖2 修改控制框架

    Fig.3 Calculate the new global coordinates圖3 計(jì)算新全局坐標(biāo)

    In the FFD control frame,the axes of local coordinate system are along the directions of length,width and height.Assuming the control points are uniformly selected along each edge,then[8-9]

    Where,X0( o')is the coordinate of the local coordinate system’s origin point in the Cartesian coordinates.X is the Cartesian coordinate of an arbitrary point inside the control frame,whose local coordinates are s,t,u separately.The value range of s,t,u is 0≤s,t,u≤1.S,T,U are the axis vectors of local coordinate system,and l+1,m+1,n+1 are the numbers of control points of the control frame along S,T,U directions.

    Since the basis function of the FFD algorithm chosen in present analysis is NURBS function,the method has several geometric characteristics similar to NURBS surface,such as continuity,convex hull property and local approximation,therefore,this method is more suitable for afterbody geometry design compared to other methods.In addition,the FFD method will guarantee the derivation of an arbitrary order’s continuity and control the changing degree of volume.

    2 Infinite interpolation deforming grid technique

    In the application of aerodynamics optimization,the deforming grid technique is required for efficiency as the shape changing.In present work,grid-point-connecting multi-block structured grid is used to discretize computational domain.The number of control frames could be set hundreds or thousands when the shape is complex.The requirements of the deforming grid technique are as follows:

    1)The capacity to express a new shape.

    2)The high quality of grid.

    3)Using parallel computing method to guarantee the efficiency.

    4)Guarantee the characteristics of grid point connection.

    In order to get a rational multi-block’s topological structure when the shape is changed,the volume spline interpolation technique is used to compute the deformation of block vertices.The technique is expressed as follows:

    After the displacement is computed by volume spline interpolation technique,other grid points’coordinates in the block could be computed by transfinite interpolation(TFI)method[10], which includes three steps of iteration and has been widely used to achieve grid deformation.The first step is by linear interpolation to compute the inner displacement along theξdirection.

    The second step is to superimpose the displacement to theηdirection

    Similarly,along theγdirection

    NI,NJ,NK are the dimensions of the grid points in the grid blocks.

    The expression of Sξ,η,γis

    Similarly the tξ,η,γ,uξ,η,γcould be computed.Thus,the grid deformation(ξ,η,γ)=(ξ,η,γ)is achieved.

    3 Flow field numerical simulation

    The flow control equations are 3-dimensional compressible unsteady Navier-Stokes equations in an integral form.The expression in Cartesian coordinate system is

    The turbulence model considered here is Menter’s k-ω SSTmodel.For spatial discretion,3-order upwind MUSCL interpolation ROE scheme is used,together with multi-grid and parallel computing technique[11]. By comparing the numerical results with experimental results of DLR-F6 wing-body as shown in Figure 4,the CFD code used in this study is proved to be reliable.The computational condition is Ma∞=0.75,Re=3.0×106

    Fig.4 DLR-F6 wing-body lift-drag performance圖4 DLR-F6升阻極曲線

    4 Optimization framework

    The modified particle swarm algorithm has been applied as the optimization method.Sun J and his fellows put forward the Quantum-Behaved Particle Swarm Optimization,namely the QPSO[12].It is different from SPSO that the searching pattern is along the track,the position of the searching particle is determined by the probability density function,which leads to a better global searching performance.The equation of the QPSO algorithm can be denoted as follows

    In the functions above,M is the number of particles of the swarm population,mbest is the average position of the pbest of each particle,φ,u are random numbers in(0,1),which are selected as the probability of both+50%and-50%,andβis the elastic coefficient as the control parameter in QPSO,with its value decreasing from 1.0 to 0.5 as the iteration continues.

    As high accuracy surrogate model is the key technique for improving the efficiency of optimization design,the surrogate model used here is the modified Kriging surrogate model.Kriging surrogate model originates from the spatial statistics in geography.It is the unbiased estimation model whose estimating variance is the smallest.Kriging surrogate model has the characteristics of local approximation by the correlation function[13],can well predict the function value distribution at the unknown points.The relation of response and design variables can be denoted as the following equation in Kriging surrogate model,

    Regressive model F(x)is the universal approximation of the design space.It is the certainty part,and can be divided into three categories:0 order(constant),1-order(linear)and 2-order(binomial).z(x)is a statistical random procedure whose average value is 0 and variance is σ2.Covariance of two interpolation points is,

    Where R is the correlation function of point x(i)and x(j).The Gaussian Function is used in present analysis and also in most applications,

    Therelated model parameterθkis determined by maximum similarity approximation,

    The related model parameter has great influence on the performance of the surrogate model.The traditional Kriging method solves the related parametersusing pattern searching method,which depends on the selection of the initial points and is easy to be trapped in the local optimum area.The standard particle swarm algorithm is used to optimize the related parameter of the Kriging surrogate model in order to improve the approximation[14].

    For the reason that uniform designcan describe the sample space characteristics commendably and uniform samples have certain advantages in solving multi-level problem,the samples in present analysis are selected uniformly.the average relative error(ARE)is compared between the initial Kriging model and the modified Kriging model by PSO algorithm.

    Table 1 shows that the average relative error of wing-body’s drag coefficient reduces 0.173123%after optimization and we can see from Figure 5 that the most optimized predicting values are closer to the CFD results than the values before optimization.

    表1 ARE對(duì)比Tables 1 ARE comparing

    Fig.5 The predictive compare of test sample圖5 測(cè)試樣本預(yù)測(cè)對(duì)比

    The loosening surrogate model management frame is employed in this optimization procedure[15].

    5 Numerical results

    In order to optimizea transport’s afterbody under wing interference in drag reduction,the designed afterbody shapes with and without wing interference are analyzed and compared.The grid is shown in Figure 6.Structured grid with 30 blocks in the total flow field has been used,and the grid number is 3 million.Parallel computing technique,Roe spatial discretion method,LU-SGSimplicit time advancing method,k-ω SST turbulence model and multi-grid accelerating technique are applied in the CFD calculation.The design status is:Ma∞=0.85,Re=1.0×107,and CLis fixed to the value equal to cruise lift coefficient.

    The optimization target is to reduce the drag by optimizing the configuration of afterbody at cruising status.Since the upswept angle would affect the tail down angle,a constraint is given for upswept angle to make sure it would not decrease.Considering the requirement of minimum space for capacity and body structure,the area of three sections are constrained to larger than specific value.

    The optimization design problem could be denoted as below:

    Tsectionconstraints of afterbody come from two places,as shown in Figure 7.One is the combination of floor and fuselage structure at the 78%length of the fuselage.The keel line of the fuselage must be lower than the plane of floor by 200 mm to save enough space for structure.The other is the cargo height,which must be higher than the initial height.S is the cross section area,and three cross section area is concerned at 70%,80%,90%length of the fuselage respectively.CMis the pitching moment constraint.All the constraints are added to the objective function by using linear penalty function.The objective function shows as follows

    ωiis the weight of the CDand other constraints.From the function(23),it can be found that the constraints term would equal to zero when the constraints are satisfied,otherwise,the objective would increase.Figure 8 shows the control frame and the afterbody geometry.The control point of I0and J0move rigidly to ensure that the surface’s curvature between the middlebody and afterbody are continuous.K0and Kmaxmove along Z direction that can modify the crown line and the keel line of the airplane,and Jmaxpoint moves along Y direction that can modify the maximum half breadth shapes.

    Fig.6 Grid圖6 表面網(wǎng)格

    Fig.7 Height constraints圖7 高度約束

    Fig.8 FFD control frame and afterbody geometry圖8 FFD控制框架及后體幾何外形

    The parametric method is arbitrary spatial FFD approach.The geometry surface of afterbody is deformed by changing the FFD control frame vertices as introduced previously,then the geometry parameters of afterbody,such as the section shape and the upswept angle,would be modified.The optimization algorithm is the quantum particle swarm algorithm with a population of 90.The Latin hyper-cube method has been implemented in sample selection for establishment of the surrogate model.The total number of samples is 300.

    With the Kriging surrogate model,the approximation error(%RMSE)of drag is 0.6%by crossed testification method.The optimization search is carried out for 60 generations,and the surrogate model is updated by CFD solver in every 10 generations.So the program calls updated geometry model for 6 times,and the CFD solver has been invoked 306 times.The total computational time is 734.4 machine-hours for computer with i7 3820 CPU and 8G RAM.

    The body section area at 60%of fuselage length before and after the optimization is compared in the Figure 9.The initial body section is a dual circle shape,and the optimization result decreases the curvature variance gradient while the curvature radius is increased.The section shape is no longer a dual circle,which decreases the circumferential pressure gradient and destabilization of afterbody flow.

    Fig.9 X/C=0.6 section shape comparison圖9 X/C=0.6截面外形優(yōu)化結(jié)果對(duì)比

    The symmetry contour of the body is shown in Figure 10 with the comparison of with and without the wing interference.Because the local angle of attack of the body has been reduced by the downwash of the wing,the highest location of the symmetry contour under the wing interference is more ahead than that without wing interference.The bottom contour is lower than the result without effect of wing to reduce the pressure recovery gradient.

    The pressure distributions of the plane of symmetry and plane at 60%of fuselage length are shown in Figure 11 and Figure 12,in which the decrease of the pressure recovery gradient after optimization can be clearly seen,and the pressure drag is reduced as a result.The cruising state aerodynamic characteristics before and after the optimization design are shown in table 2,in which the drag coefficient decreases by 6 counts and the liftdrag ratio increases by 3%,which is a great benefit in drag reduction.Cmis the pitching moment coefficient,which is slightly decreased to meet the constraint.From the limiting stream line on the afterbody in figure 13,we can see that the separation reduces to a small area,which would vanish if consider an APU system.

    Fig.10 Afterbody symmetry contour圖10 后體輪廓線

    Fig.11 Pressure distribution on the symmetry plane圖11 對(duì)稱面壓力分布

    Fig.12 Z/C=0.6 pressure distribution圖12 Z/C=0.6壓力分布

    表2 優(yōu)化前后氣動(dòng)特性Table 2 Aerodynamic characteristics after the optimization

    Fig.13 Limiting stream line on afterbody圖13 表面極限流線對(duì)比

    6 Conclusion

    The presented work in this paper is the optimization design of transport afterbody by FFD parametric method integrated with Kriging surrogate model and quantum particle swarm algorithm,which are adopted to establish the aerodynamic optimization design management frame.

    A numerical test is done for a typical transport afterbody considering the wing wake flow,from which the result shows that the downwash of the wing could lead to a decrease of the local angle of attack of the body flow field,so the maximum height location of the upper outline of the optimized configuration is more ahead than that without the wing interference,and the bottom outline with the wing interference becomes lower than that without the wing interference.This test shows that the afterbody flow field is affected by the wing wake flow dramatically,so the design of the afterbody should take that influence into consideration.

    The body pressure dragis reduced after the optimization design by reducing the pressure recovery gradient.The drag coefficient decreases by 6 counts with comparison of that prior optimization and the lift-drag ratio increases by 3%.

    The aerodynamic optimization design system for transport’s afterbody established in this study is of good optimization design efficiency and indicates a promising future of engineering application.

    [1]Wortman A.Reduction of fuselage form drag by vortex flows[J].Journal of Aircraft,1999,36(3):501-506.

    [2]Epatein R J,Carbonaro M C,Caudrom F.An experimental investigation of the flow field about an upswept afterbody[J].Journal of Aircraft,1994,31(6):1281-1290.

    [3]Huang Y,Ghia U,Osswald G A,et al.Analysis and numerical simulation of 3-D flow past axisymmetric afterbody using Navier-Stokes equations[R].AIAA 1993-0683.

    [4]Zhang Binqian,Wang Yuanyuan Duan Zhuoyi,et al.Design method for large upswept afterbody of transport aircraft[J].Acta Aeronautica et Astronautica Sinica,2010,31(10):1933-1939.(in Chinese)張彬乾,王元元,段卓毅,等.大上翹機(jī)身后體設(shè)計(jì)方法[J].航空學(xué)報(bào),2010,31(10):1933-1939.

    [5]Kong Fanmei,Hua Jun.Effects of geometry parameters and flow parameters on drag coefficient of upswept afterbodies[J].Journal of Beijing University of Aeronautics and Astronautics,2003,29(1):39-42.

    [6]Kolesar C E,May F.An after drag prediction technique for military airlifters[R].AIAA 1983-1787.

    [7]Thomas J Otahal,et al.An investigation of two dimensional CAD generated models with body decoupled cartesian grids for DSMC[C]//34thAIAA Thermophysics Conference,2000:19-22.

    [8]Zhu Xinxiong.The technical of free curve and surface sculpt[M].Beijing:Science Press,2000.(in Chinese)朱雄心.自由曲線曲面造型技術(shù)[M].北京:科學(xué)出版社,2000.

    [9]Sederberg T W,Parry SR.Freeform deformation of solid geometric models[J].Computer Graphics,19886,22(4):151-160.

    [10]Smith R E.Transfinite interpolation(TFI)generation systems[M].eds.N.P.Weatherill,J.F.Thompson,B.K.Soni.Handbook of Grid Generation,CRC Press,1999.

    [11]Menter F R.Two-equation eddy-viscosity turbulence models for engineering applications[J].AIAA Journal,1994,32(8):269-289.

    [12]Chen Peng,Li Jian,Guan Tao.The optimization of parameters of kriging correlation model based on particle swarm optimization[J].Microelectronics& Computer,2009,26(4):178-181.(in Chinese)陳鵬,李劍,管濤.基于PSO的Kriging相關(guān)模型參數(shù)優(yōu)化[J].微電子學(xué)與計(jì)算機(jī),2009,26(4):178-181.

    [13]Shinkyu Jeong,Mitsuhiro Murayama,Kazuomi Yamamoto.Efficient optimization design method using Kriging model[J].Journal of AIAA,2005,42(2):413-420.

    [14]Sun J,F(xiàn)eng B,Xu W B.Particle swarm optimization with particles having quantum behavior[C]//Proc of the IEEE Congress on Evolutionary Computation,2004:325-331.

    [15]Knill D L,Giunta A A.Response surface models combining linear and euler aerodynamics for supersonic transport design[J].Journal of Aircraft,1999,36(1):75-86.

    猜你喜歡
    張彬外形框架
    Compared discharge characteristics and film modifications of atmospheric pressure plasma jets with two different electrode geometries
    EXISTENCE RESULTS FOR SINGULAR FRACTIONAL p-KIRCHHOFF PROBLEMS*
    比外形,都不同
    框架
    廣義框架的不相交性
    酒鬼報(bào)警
    WTO框架下
    法大研究生(2017年1期)2017-04-10 08:55:06
    張彬 作品選登
    論袁牧之“外形的演技”
    一種基于OpenStack的云應(yīng)用開發(fā)框架
    国产视频一区二区在线看| 精品国产超薄肉色丝袜足j| 久久久久亚洲av毛片大全| 日本 av在线| 国产一区二区三区在线臀色熟女| 69精品国产乱码久久久| 日韩大尺度精品在线看网址 | 久热这里只有精品99| 一a级毛片在线观看| 长腿黑丝高跟| 国产亚洲av嫩草精品影院| 成人18禁在线播放| 亚洲成人免费电影在线观看| 欧美老熟妇乱子伦牲交| 国产亚洲av嫩草精品影院| 级片在线观看| 久久精品91蜜桃| 两个人免费观看高清视频| 狠狠狠狠99中文字幕| 国产av又大| 非洲黑人性xxxx精品又粗又长| 丰满人妻熟妇乱又伦精品不卡| 色婷婷久久久亚洲欧美| 久热这里只有精品99| 亚洲一码二码三码区别大吗| 亚洲国产精品久久男人天堂| 欧美日韩瑟瑟在线播放| 日韩欧美一区视频在线观看| 日韩视频一区二区在线观看| 免费看十八禁软件| 久久影院123| 变态另类成人亚洲欧美熟女 | 欧美久久黑人一区二区| 免费在线观看影片大全网站| 天堂√8在线中文| 久久九九热精品免费| 成人18禁高潮啪啪吃奶动态图| 一级黄色大片毛片| 91麻豆精品激情在线观看国产| 精品免费久久久久久久清纯| 少妇的丰满在线观看| 色综合婷婷激情| 精品电影一区二区在线| 俄罗斯特黄特色一大片| 一卡2卡三卡四卡精品乱码亚洲| 免费在线观看视频国产中文字幕亚洲| 在线观看免费视频网站a站| 满18在线观看网站| www.精华液| 激情在线观看视频在线高清| 在线永久观看黄色视频| 午夜精品在线福利| 精品高清国产在线一区| 涩涩av久久男人的天堂| 欧美日韩乱码在线| 成人18禁在线播放| 日韩精品中文字幕看吧| 少妇粗大呻吟视频| videosex国产| 18禁裸乳无遮挡免费网站照片 | 色播在线永久视频| 精品国内亚洲2022精品成人| 精品卡一卡二卡四卡免费| 午夜日韩欧美国产| 伦理电影免费视频| 久久午夜综合久久蜜桃| 国产激情欧美一区二区| 欧美日本视频| 亚洲av五月六月丁香网| 国产午夜福利久久久久久| 夜夜看夜夜爽夜夜摸| 亚洲 国产 在线| 国产亚洲精品第一综合不卡| 午夜福利一区二区在线看| 90打野战视频偷拍视频| 夜夜爽天天搞| 国产xxxxx性猛交| av在线天堂中文字幕| 桃色一区二区三区在线观看| 97人妻天天添夜夜摸| 国产精品久久久久久精品电影 | 亚洲免费av在线视频| 好看av亚洲va欧美ⅴa在| 色在线成人网| 最近最新中文字幕大全免费视频| 久久久久精品国产欧美久久久| 久热爱精品视频在线9| 一区二区日韩欧美中文字幕| 中文字幕另类日韩欧美亚洲嫩草| 操出白浆在线播放| 欧美精品啪啪一区二区三区| 亚洲av美国av| 日本一区二区免费在线视频| 超碰成人久久| 日韩 欧美 亚洲 中文字幕| 久久久久久大精品| 黄片大片在线免费观看| 久久久久久久精品吃奶| 在线观看午夜福利视频| 久久性视频一级片| 女人被狂操c到高潮| 中文字幕色久视频| 亚洲avbb在线观看| 亚洲国产精品久久男人天堂| 高潮久久久久久久久久久不卡| 亚洲国产日韩欧美精品在线观看 | 亚洲精品国产色婷婷电影| 欧美最黄视频在线播放免费| 在线观看舔阴道视频| 午夜两性在线视频| 色在线成人网| 露出奶头的视频| 亚洲一区高清亚洲精品| 欧美日韩精品网址| 国产三级黄色录像| 欧美成狂野欧美在线观看| 午夜久久久久精精品| 欧美日韩精品网址| 国产亚洲欧美在线一区二区| 91麻豆av在线| 亚洲人成网站在线播放欧美日韩| or卡值多少钱| 校园春色视频在线观看| 国产伦一二天堂av在线观看| 亚洲精品中文字幕在线视频| 日韩欧美免费精品| 天堂影院成人在线观看| 99国产极品粉嫩在线观看| а√天堂www在线а√下载| 国产成人啪精品午夜网站| 亚洲中文字幕日韩| 亚洲性夜色夜夜综合| 自拍欧美九色日韩亚洲蝌蚪91| 九色亚洲精品在线播放| 国产精品免费一区二区三区在线| 国产精品久久久久久人妻精品电影| 中文字幕人成人乱码亚洲影| 欧美成人一区二区免费高清观看 | 热re99久久国产66热| 亚洲国产高清在线一区二区三 | 亚洲专区国产一区二区| 美女免费视频网站| 深夜精品福利| 国产精品九九99| 日日干狠狠操夜夜爽| 国产成人av教育| 最近最新中文字幕大全电影3 | 麻豆一二三区av精品| 超碰成人久久| 国产一级毛片七仙女欲春2 | 90打野战视频偷拍视频| 桃色一区二区三区在线观看| 成人特级黄色片久久久久久久| av有码第一页| 级片在线观看| 91成人精品电影| 日本 欧美在线| 淫妇啪啪啪对白视频| 久久欧美精品欧美久久欧美| 丝袜人妻中文字幕| 色综合婷婷激情| av在线播放免费不卡| 禁无遮挡网站| 欧美日本视频| 色综合亚洲欧美另类图片| 欧美日韩一级在线毛片| 777久久人妻少妇嫩草av网站| 91精品三级在线观看| 国产区一区二久久| 久久婷婷成人综合色麻豆| 丝袜美腿诱惑在线| 91大片在线观看| 伊人久久大香线蕉亚洲五| 欧美成人性av电影在线观看| 亚洲在线自拍视频| 国产野战对白在线观看| 给我免费播放毛片高清在线观看| 欧美国产精品va在线观看不卡| 欧美老熟妇乱子伦牲交| 国产97色在线日韩免费| 国产精品国产高清国产av| 国产成年人精品一区二区| 国产麻豆成人av免费视频| 欧美国产日韩亚洲一区| 国产精品 国内视频| 极品人妻少妇av视频| netflix在线观看网站| 色老头精品视频在线观看| 天天躁狠狠躁夜夜躁狠狠躁| 中文字幕另类日韩欧美亚洲嫩草| 欧美人与性动交α欧美精品济南到| 最近最新中文字幕大全电影3 | 国产精品秋霞免费鲁丝片| 看黄色毛片网站| 午夜影院日韩av| 国产午夜精品久久久久久| 非洲黑人性xxxx精品又粗又长| 久久精品亚洲精品国产色婷小说| 免费看a级黄色片| 黑丝袜美女国产一区| aaaaa片日本免费| 久久精品国产综合久久久| 每晚都被弄得嗷嗷叫到高潮| 久久久国产精品麻豆| 97超级碰碰碰精品色视频在线观看| 搡老熟女国产l中国老女人| 99国产极品粉嫩在线观看| 国产精品av久久久久免费| 国产成人欧美在线观看| 欧美精品啪啪一区二区三区| 天天躁狠狠躁夜夜躁狠狠躁| 国产精品日韩av在线免费观看 | 黑人巨大精品欧美一区二区mp4| 日韩欧美三级三区| 午夜亚洲福利在线播放| 国产国语露脸激情在线看| 国产高清videossex| 色播在线永久视频| 精品电影一区二区在线| 99在线人妻在线中文字幕| 久久久久国内视频| 18禁黄网站禁片午夜丰满| 人人妻人人澡人人看| 午夜久久久在线观看| 黄色a级毛片大全视频| 久久婷婷成人综合色麻豆| 99热只有精品国产| 久9热在线精品视频| 日韩欧美国产在线观看| 老司机深夜福利视频在线观看| 一区二区三区激情视频| 日本 av在线| 18美女黄网站色大片免费观看| 午夜免费鲁丝| 久久亚洲精品不卡| 中文字幕色久视频| 免费看a级黄色片| 两性夫妻黄色片| 日韩精品中文字幕看吧| 免费看美女性在线毛片视频| 人人澡人人妻人| 搡老岳熟女国产| 欧美日韩福利视频一区二区| 亚洲电影在线观看av| 性色av乱码一区二区三区2| 国内久久婷婷六月综合欲色啪| 日本五十路高清| av在线天堂中文字幕| 久久久久久久久免费视频了| 日韩欧美一区二区三区在线观看| 99久久综合精品五月天人人| 欧美日本亚洲视频在线播放| 欧美人与性动交α欧美精品济南到| 欧美日韩精品网址| 两个人视频免费观看高清| 欧美另类亚洲清纯唯美| 美女免费视频网站| 在线av久久热| 女人高潮潮喷娇喘18禁视频| 日韩有码中文字幕| 亚洲一区二区三区色噜噜| 国产成人精品在线电影| 黄色a级毛片大全视频| 亚洲精品久久成人aⅴ小说| 97超级碰碰碰精品色视频在线观看| 999久久久精品免费观看国产| 男女下面插进去视频免费观看| 亚洲中文日韩欧美视频| 天堂√8在线中文| 欧美国产日韩亚洲一区| 亚洲欧美日韩无卡精品| 最新在线观看一区二区三区| 国产高清视频在线播放一区| 真人一进一出gif抽搐免费| 美女国产高潮福利片在线看| 国产成人精品在线电影| 真人做人爱边吃奶动态| 亚洲精品国产精品久久久不卡| 正在播放国产对白刺激| 国产免费av片在线观看野外av| 国产区一区二久久| 精品人妻1区二区| 欧美成狂野欧美在线观看| 九色亚洲精品在线播放| 啪啪无遮挡十八禁网站| 亚洲午夜精品一区,二区,三区| 91成年电影在线观看| 精品国产乱子伦一区二区三区| 在线视频色国产色| 99久久国产精品久久久| 大型黄色视频在线免费观看| 真人一进一出gif抽搐免费| 变态另类成人亚洲欧美熟女 | 日韩中文字幕欧美一区二区| 日本 欧美在线| 中文字幕精品免费在线观看视频| 亚洲人成77777在线视频| 午夜免费鲁丝| 十八禁网站免费在线| www.熟女人妻精品国产| 亚洲成人精品中文字幕电影| 精品久久久久久久毛片微露脸| 亚洲成国产人片在线观看| 国产乱人伦免费视频| 国内毛片毛片毛片毛片毛片| 91在线观看av| 成在线人永久免费视频| 咕卡用的链子| 极品教师在线免费播放| 欧美日本视频| 国产精品av久久久久免费| 久久国产亚洲av麻豆专区| 亚洲成av片中文字幕在线观看| 国内久久婷婷六月综合欲色啪| 免费少妇av软件| 99久久精品国产亚洲精品| 国产精品综合久久久久久久免费 | 成人国产综合亚洲| cao死你这个sao货| 曰老女人黄片| 91av网站免费观看| 老司机深夜福利视频在线观看| 日韩视频一区二区在线观看| 在线观看舔阴道视频| 成人亚洲精品一区在线观看| 国产成人av激情在线播放| 18禁国产床啪视频网站| 久久久久久亚洲精品国产蜜桃av| 女人高潮潮喷娇喘18禁视频| 最新在线观看一区二区三区| 欧美av亚洲av综合av国产av| 亚洲欧美精品综合久久99| 亚洲熟女毛片儿| 国产一区二区三区在线臀色熟女| 午夜日韩欧美国产| 黄色视频,在线免费观看| 色综合婷婷激情| 国产精品日韩av在线免费观看 | 午夜亚洲福利在线播放| 不卡av一区二区三区| 91av网站免费观看| 国产精品乱码一区二三区的特点 | 国产精品99久久99久久久不卡| 99国产综合亚洲精品| 男男h啪啪无遮挡| 亚洲精品国产精品久久久不卡| 看免费av毛片| 婷婷六月久久综合丁香| 热99re8久久精品国产| 亚洲成人久久性| 成年人黄色毛片网站| 亚洲成人免费电影在线观看| 制服诱惑二区| 老熟妇仑乱视频hdxx| 久9热在线精品视频| 少妇 在线观看| 国产伦人伦偷精品视频| 中文字幕高清在线视频| 一个人免费在线观看的高清视频| 精品久久久精品久久久| 女性生殖器流出的白浆| 成人18禁高潮啪啪吃奶动态图| 美女大奶头视频| 免费在线观看黄色视频的| 97人妻天天添夜夜摸| 欧美日韩亚洲综合一区二区三区_| netflix在线观看网站| 免费高清在线观看日韩| 成人永久免费在线观看视频| 日韩视频一区二区在线观看| 最近最新中文字幕大全电影3 | www.999成人在线观看| 色尼玛亚洲综合影院| 日韩免费av在线播放| 中文字幕人妻熟女乱码| 成人国产一区最新在线观看| 99久久国产精品久久久| 午夜影院日韩av| www国产在线视频色| 久久欧美精品欧美久久欧美| 亚洲熟妇熟女久久| 午夜日韩欧美国产| 国产亚洲精品久久久久5区| 欧美乱色亚洲激情| 国产成人精品无人区| 亚洲午夜精品一区,二区,三区| 日韩av在线大香蕉| 侵犯人妻中文字幕一二三四区| 1024视频免费在线观看| 亚洲欧美激情在线| 如日韩欧美国产精品一区二区三区| 999精品在线视频| 亚洲一卡2卡3卡4卡5卡精品中文| 亚洲全国av大片| 黄色片一级片一级黄色片| 日韩精品青青久久久久久| 日韩一卡2卡3卡4卡2021年| 国产三级在线视频| 成熟少妇高潮喷水视频| 欧美午夜高清在线| 91老司机精品| 欧美激情极品国产一区二区三区| 韩国av一区二区三区四区| 91av网站免费观看| 欧美日韩精品网址| 成人手机av| 99久久久亚洲精品蜜臀av| 国产国语露脸激情在线看| av天堂在线播放| 国产国语露脸激情在线看| 午夜免费成人在线视频| 国产精品亚洲av一区麻豆| 丁香欧美五月| 久久伊人香网站| 国产精品久久久久久精品电影 | 国产单亲对白刺激| 在线观看www视频免费| 极品教师在线免费播放| 啦啦啦观看免费观看视频高清 | 久久久国产成人精品二区| 麻豆久久精品国产亚洲av| 1024香蕉在线观看| 久久这里只有精品19| 麻豆久久精品国产亚洲av| 十八禁人妻一区二区| 波多野结衣巨乳人妻| 一a级毛片在线观看| 色av中文字幕| 最好的美女福利视频网| 日本撒尿小便嘘嘘汇集6| 男人舔女人的私密视频| 欧美色视频一区免费| 在线十欧美十亚洲十日本专区| 国产成人精品在线电影| 欧美一级a爱片免费观看看 | 久久久久久久久免费视频了| 亚洲国产毛片av蜜桃av| 99精品久久久久人妻精品| 国产精品98久久久久久宅男小说| 久久婷婷人人爽人人干人人爱 | 岛国在线观看网站| 日本精品一区二区三区蜜桃| 午夜福利成人在线免费观看| 成人av一区二区三区在线看| 大型av网站在线播放| 国产免费男女视频| 天堂动漫精品| 丰满人妻熟妇乱又伦精品不卡| 国产精品影院久久| 麻豆国产av国片精品| 成年版毛片免费区| 男女之事视频高清在线观看| 国产激情欧美一区二区| 亚洲激情在线av| 国内久久婷婷六月综合欲色啪| av欧美777| 成人三级做爰电影| 一级a爱片免费观看的视频| 亚洲男人的天堂狠狠| 欧美丝袜亚洲另类 | 国产精品 欧美亚洲| 欧美黄色片欧美黄色片| 女警被强在线播放| 国产xxxxx性猛交| 国产片内射在线| 亚洲av美国av| 91九色精品人成在线观看| 久久久久久久精品吃奶| 不卡一级毛片| 最新在线观看一区二区三区| 在线永久观看黄色视频| 久久国产亚洲av麻豆专区| 性色av乱码一区二区三区2| 麻豆成人av在线观看| 别揉我奶头~嗯~啊~动态视频| 亚洲av成人不卡在线观看播放网| 自拍欧美九色日韩亚洲蝌蚪91| 国产精品,欧美在线| 日日夜夜操网爽| 国产伦人伦偷精品视频| 欧美成人午夜精品| 很黄的视频免费| av欧美777| 日韩大尺度精品在线看网址 | 777久久人妻少妇嫩草av网站| 久久午夜亚洲精品久久| 国产三级在线视频| 国产人伦9x9x在线观看| 90打野战视频偷拍视频| 午夜激情av网站| 麻豆国产av国片精品| 亚洲成人精品中文字幕电影| 久久午夜综合久久蜜桃| 久久久久久大精品| 国产精品99久久99久久久不卡| 给我免费播放毛片高清在线观看| 免费看a级黄色片| 国内久久婷婷六月综合欲色啪| 在线天堂中文资源库| 午夜福利影视在线免费观看| 日日干狠狠操夜夜爽| 午夜激情av网站| 中国美女看黄片| 国产亚洲精品久久久久久毛片| 国内毛片毛片毛片毛片毛片| 18美女黄网站色大片免费观看| 搞女人的毛片| 叶爱在线成人免费视频播放| 麻豆一二三区av精品| 亚洲av成人不卡在线观看播放网| 长腿黑丝高跟| 伦理电影免费视频| 婷婷六月久久综合丁香| 欧美日韩亚洲国产一区二区在线观看| 丁香欧美五月| 精品午夜福利视频在线观看一区| 巨乳人妻的诱惑在线观看| 两人在一起打扑克的视频| 人人妻人人澡欧美一区二区 | 精品久久久精品久久久| 精品国产乱码久久久久久男人| 91九色精品人成在线观看| 免费高清在线观看日韩| av电影中文网址| 亚洲七黄色美女视频| 亚洲一区中文字幕在线| 色播亚洲综合网| 亚洲五月色婷婷综合| 久久性视频一级片| 丝袜在线中文字幕| 操出白浆在线播放| 精品人妻1区二区| 黄片大片在线免费观看| 亚洲人成网站在线播放欧美日韩| 黄片大片在线免费观看| 麻豆久久精品国产亚洲av| 欧美人与性动交α欧美精品济南到| 国产精品综合久久久久久久免费 | 精品午夜福利视频在线观看一区| 纯流量卡能插随身wifi吗| 国产麻豆69| 伊人久久大香线蕉亚洲五| 桃色一区二区三区在线观看| 国产三级在线视频| 性欧美人与动物交配| 美女午夜性视频免费| 国产精品 欧美亚洲| videosex国产| 在线观看免费视频日本深夜| 免费不卡黄色视频| 18禁国产床啪视频网站| 无限看片的www在线观看| 久久亚洲真实| 此物有八面人人有两片| 日日摸夜夜添夜夜添小说| 97人妻天天添夜夜摸| 国产一区二区三区综合在线观看| 精品国内亚洲2022精品成人| 久久精品人人爽人人爽视色| 好看av亚洲va欧美ⅴa在| 国产伦人伦偷精品视频| 久久久久国产精品人妻aⅴ院| 亚洲aⅴ乱码一区二区在线播放 | 久久久久国内视频| 亚洲精品国产区一区二| 人人妻人人澡人人看| 巨乳人妻的诱惑在线观看| 国产人伦9x9x在线观看| 国产欧美日韩精品亚洲av| 精品高清国产在线一区| 91av网站免费观看| 国产乱人伦免费视频| 午夜福利视频1000在线观看 | 18禁黄网站禁片午夜丰满| 涩涩av久久男人的天堂| 久久久久亚洲av毛片大全| 亚洲免费av在线视频| av在线播放免费不卡| 男男h啪啪无遮挡| 99riav亚洲国产免费| 国产在线观看jvid| 久久国产亚洲av麻豆专区| 在线观看午夜福利视频| 亚洲国产精品sss在线观看| 日日夜夜操网爽| 久久天躁狠狠躁夜夜2o2o| 人人妻人人澡欧美一区二区 | 欧美激情高清一区二区三区| 久久久久久免费高清国产稀缺| 亚洲av电影不卡..在线观看| 久久精品人人爽人人爽视色| 一区二区三区激情视频| 午夜福利成人在线免费观看| 亚洲人成77777在线视频| 妹子高潮喷水视频| www.精华液| 叶爱在线成人免费视频播放| 悠悠久久av| 成人18禁在线播放| 国产伦人伦偷精品视频| 18禁观看日本| 精品第一国产精品| 亚洲一码二码三码区别大吗| 神马国产精品三级电影在线观看 | 日本三级黄在线观看| 久久久久久久午夜电影| 亚洲中文字幕日韩| 成人亚洲精品av一区二区| 嫩草影院精品99| 黄色 视频免费看| 久久久久精品国产欧美久久久| 搡老妇女老女人老熟妇| 在线观看www视频免费| 日本vs欧美在线观看视频| 在线永久观看黄色视频|