• <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ā)框架
    大型av网站在线播放| 成人午夜高清在线视频 | 给我免费播放毛片高清在线观看| 黑丝袜美女国产一区| 亚洲性夜色夜夜综合| 变态另类丝袜制服| 国产av不卡久久| 18禁国产床啪视频网站| 国产av一区二区精品久久| 国产欧美日韩一区二区三| 久久精品国产清高在天天线| 久久国产精品男人的天堂亚洲| 2021天堂中文幕一二区在线观 | 青草久久国产| 俺也久久电影网| 久久中文字幕人妻熟女| 国产亚洲精品久久久久5区| 久久精品国产亚洲av高清一级| 欧美激情高清一区二区三区| 大型av网站在线播放| 18禁裸乳无遮挡免费网站照片 | 最新在线观看一区二区三区| 夜夜夜夜夜久久久久| 桃色一区二区三区在线观看| 亚洲第一欧美日韩一区二区三区| 午夜福利成人在线免费观看| svipshipincom国产片| 亚洲成国产人片在线观看| 夜夜夜夜夜久久久久| 十分钟在线观看高清视频www| 亚洲免费av在线视频| 午夜福利成人在线免费观看| 精品日产1卡2卡| 久久人妻av系列| 精品久久久久久久人妻蜜臀av| 久久精品亚洲精品国产色婷小说| 国产精品亚洲美女久久久| av免费在线观看网站| 激情在线观看视频在线高清| 国产精品98久久久久久宅男小说| 亚洲成人免费电影在线观看| 国产一区二区三区视频了| 一进一出抽搐gif免费好疼| 亚洲精品中文字幕在线视频| 亚洲成人久久性| 久久久久久人人人人人| 99国产极品粉嫩在线观看| 婷婷精品国产亚洲av在线| 精品人妻1区二区| 欧美最黄视频在线播放免费| 久久九九热精品免费| 成年人黄色毛片网站| 91国产中文字幕| 无人区码免费观看不卡| 麻豆国产av国片精品| 欧美日韩一级在线毛片| 亚洲成人久久爱视频| 非洲黑人性xxxx精品又粗又长| aaaaa片日本免费| 天天添夜夜摸| 一个人免费在线观看的高清视频| 黑丝袜美女国产一区| 免费观看人在逋| cao死你这个sao货| 一区二区三区精品91| 国产极品粉嫩免费观看在线| 女人被狂操c到高潮| 狂野欧美激情性xxxx| 亚洲精品国产一区二区精华液| 美女扒开内裤让男人捅视频| 欧美中文综合在线视频| 高潮久久久久久久久久久不卡| 亚洲精品av麻豆狂野| 欧美乱色亚洲激情| 手机成人av网站| 欧美成人午夜精品| 日韩免费av在线播放| 搡老岳熟女国产| 久久天堂一区二区三区四区| 哪里可以看免费的av片| 99国产综合亚洲精品| 久久精品国产亚洲av高清一级| 99国产精品一区二区蜜桃av| 18禁裸乳无遮挡免费网站照片 | 久久久久久久久免费视频了| 久久久久久久久中文| 色婷婷久久久亚洲欧美| 亚洲中文字幕一区二区三区有码在线看 | 手机成人av网站| 欧美日韩黄片免| 久久国产乱子伦精品免费另类| 欧美日韩中文字幕国产精品一区二区三区| 亚洲七黄色美女视频| 欧美色视频一区免费| 免费看日本二区| 午夜久久久久精精品| 久久精品国产综合久久久| 一本综合久久免费| 亚洲中文字幕一区二区三区有码在线看 | 90打野战视频偷拍视频| av超薄肉色丝袜交足视频| 日韩有码中文字幕| 国产97色在线日韩免费| 两个人看的免费小视频| av电影中文网址| 国产亚洲精品久久久久久毛片| 国产亚洲精品第一综合不卡| 长腿黑丝高跟| 亚洲七黄色美女视频| 成人手机av| 97超级碰碰碰精品色视频在线观看| 国产av一区在线观看免费| 精品久久久久久久人妻蜜臀av| 国产精品1区2区在线观看.| 99国产极品粉嫩在线观看| 88av欧美| 69av精品久久久久久| 午夜福利视频1000在线观看| www国产在线视频色| 国产精品一区二区免费欧美| 欧美乱码精品一区二区三区| 亚洲五月色婷婷综合| 亚洲自拍偷在线| 欧美激情极品国产一区二区三区| 久久精品91无色码中文字幕| 在线av久久热| 中文字幕最新亚洲高清| 女性被躁到高潮视频| 午夜福利在线在线| 88av欧美| 满18在线观看网站| 午夜精品在线福利| 在线播放国产精品三级| 精品国产国语对白av| 在线观看www视频免费| 国产成人系列免费观看| 长腿黑丝高跟| 在线观看免费日韩欧美大片| 欧美另类亚洲清纯唯美| 脱女人内裤的视频| 国产精品免费视频内射| 色在线成人网| 亚洲中文av在线| 一个人观看的视频www高清免费观看 | 人人妻人人澡人人看| 午夜久久久久精精品| av电影中文网址| 国产一级毛片七仙女欲春2 | 午夜精品久久久久久毛片777| 国产精品久久电影中文字幕| 亚洲av成人一区二区三| 亚洲成人免费电影在线观看| 成人三级黄色视频| 久久久久九九精品影院| 午夜老司机福利片| 精品第一国产精品| 中文字幕精品免费在线观看视频| 天天躁夜夜躁狠狠躁躁| 精品乱码久久久久久99久播| 一夜夜www| 18禁黄网站禁片午夜丰满| 国产av一区在线观看免费| 黄色毛片三级朝国网站| 午夜两性在线视频| 此物有八面人人有两片| 日韩欧美一区视频在线观看| 色综合站精品国产| 国产成人精品无人区| 精品一区二区三区四区五区乱码| 国产精品 欧美亚洲| 欧美日本视频| 动漫黄色视频在线观看| 在线观看一区二区三区| 免费看美女性在线毛片视频| 久久午夜亚洲精品久久| 精品日产1卡2卡| 99久久99久久久精品蜜桃| 麻豆成人av在线观看| 国产在线观看jvid| 91成年电影在线观看| 在线观看免费日韩欧美大片| 熟女少妇亚洲综合色aaa.| 两性午夜刺激爽爽歪歪视频在线观看 | 国产高清videossex| 日韩精品中文字幕看吧| 日韩欧美国产一区二区入口| 精品乱码久久久久久99久播| 国产成人精品无人区| 亚洲免费av在线视频| 久久国产亚洲av麻豆专区| 国产精品免费视频内射| 成人三级黄色视频| 午夜福利免费观看在线| 亚洲一卡2卡3卡4卡5卡精品中文| 在线播放国产精品三级| 欧美激情久久久久久爽电影| 国产在线观看jvid| 免费看美女性在线毛片视频| 日本一区二区免费在线视频| 黑人欧美特级aaaaaa片| 亚洲人成网站在线播放欧美日韩| 一区福利在线观看| 亚洲av第一区精品v没综合| 国内精品久久久久久久电影| 夜夜躁狠狠躁天天躁| 一级a爱视频在线免费观看| 无遮挡黄片免费观看| 女性生殖器流出的白浆| 欧美成狂野欧美在线观看| 最新美女视频免费是黄的| 午夜激情av网站| 99re在线观看精品视频| 久久久久久大精品| 桃红色精品国产亚洲av| 正在播放国产对白刺激| 国产精品一区二区精品视频观看| 国产又色又爽无遮挡免费看| 国产精品精品国产色婷婷| 国产极品粉嫩免费观看在线| 国产激情偷乱视频一区二区| 夜夜看夜夜爽夜夜摸| 在线观看一区二区三区| 俄罗斯特黄特色一大片| 国产日本99.免费观看| 香蕉久久夜色| 人人妻,人人澡人人爽秒播| 成人亚洲精品一区在线观看| 免费在线观看黄色视频的| 九色国产91popny在线| 中文亚洲av片在线观看爽| 香蕉av资源在线| 一区二区三区精品91| 成熟少妇高潮喷水视频| 亚洲自偷自拍图片 自拍| 日韩精品青青久久久久久| 欧美日韩中文字幕国产精品一区二区三区| 一二三四社区在线视频社区8| 青草久久国产| 国产日本99.免费观看| 级片在线观看| 国产精品1区2区在线观看.| 精品国产一区二区三区四区第35| а√天堂www在线а√下载| 久久久久久亚洲精品国产蜜桃av| 亚洲最大成人中文| 亚洲片人在线观看| 最近最新免费中文字幕在线| 欧美中文日本在线观看视频| 欧美日本视频| 999久久久国产精品视频| 欧美日韩中文字幕国产精品一区二区三区| 国产麻豆成人av免费视频| 99国产精品99久久久久| www.自偷自拍.com| 天天添夜夜摸| 国产精品亚洲av一区麻豆| 欧美精品啪啪一区二区三区| 国产在线观看jvid| cao死你这个sao货| 精品欧美一区二区三区在线| 国产精品乱码一区二三区的特点| 18禁黄网站禁片免费观看直播| 久久久精品欧美日韩精品| 日韩一卡2卡3卡4卡2021年| 亚洲在线自拍视频| 亚洲国产高清在线一区二区三 | 18美女黄网站色大片免费观看| 成人18禁在线播放| 婷婷精品国产亚洲av| 日韩精品免费视频一区二区三区| а√天堂www在线а√下载| 99久久精品国产亚洲精品| 精品久久久久久久人妻蜜臀av| 婷婷精品国产亚洲av在线| 叶爱在线成人免费视频播放| 中文亚洲av片在线观看爽| 国产男靠女视频免费网站| 丝袜美腿诱惑在线| 久久久国产成人精品二区| 婷婷精品国产亚洲av| 黄色女人牲交| 亚洲成av人片免费观看| 一级毛片高清免费大全| 99re在线观看精品视频| 这个男人来自地球电影免费观看| 中文字幕久久专区| ponron亚洲| 久久精品影院6| 黄片大片在线免费观看| 欧洲精品卡2卡3卡4卡5卡区| 国产精品亚洲美女久久久| 美女大奶头视频| 免费女性裸体啪啪无遮挡网站| 亚洲国产中文字幕在线视频| 免费无遮挡裸体视频| 成人三级黄色视频| 亚洲精品久久国产高清桃花| 欧美大码av| 亚洲中文字幕日韩| 亚洲第一青青草原| 亚洲一卡2卡3卡4卡5卡精品中文| 午夜日韩欧美国产| 精华霜和精华液先用哪个| 最近最新中文字幕大全免费视频| 亚洲成国产人片在线观看| 亚洲va日本ⅴa欧美va伊人久久| 久久精品aⅴ一区二区三区四区| 1024视频免费在线观看| 精品久久久久久久久久免费视频| 国内少妇人妻偷人精品xxx网站 | 美国免费a级毛片| 男女下面进入的视频免费午夜 | 法律面前人人平等表现在哪些方面| 国产精品1区2区在线观看.| 一二三四社区在线视频社区8| 亚洲成人久久爱视频| 中文字幕精品免费在线观看视频| 欧美日本亚洲视频在线播放| 久久中文字幕人妻熟女| 国产亚洲欧美精品永久| 亚洲熟女毛片儿| 国产精品二区激情视频| 亚洲欧美激情综合另类| 婷婷亚洲欧美| 男女床上黄色一级片免费看| 亚洲欧美激情综合另类| 欧美国产精品va在线观看不卡| 久久精品国产亚洲av香蕉五月| 高清毛片免费观看视频网站| 亚洲 欧美一区二区三区| 欧美+亚洲+日韩+国产| 淫秽高清视频在线观看| 性欧美人与动物交配| 手机成人av网站| 国产欧美日韩精品亚洲av| 人人妻人人澡人人看| 欧美乱色亚洲激情| 男女做爰动态图高潮gif福利片| 国产色视频综合| 黑人操中国人逼视频| 国产一区在线观看成人免费| 这个男人来自地球电影免费观看| 免费在线观看日本一区| 国产高清视频在线播放一区| 精品免费久久久久久久清纯| 法律面前人人平等表现在哪些方面| 精品国产亚洲在线| 国产精品爽爽va在线观看网站 | 色婷婷久久久亚洲欧美| 午夜免费观看网址| 中文字幕精品亚洲无线码一区 | 精品国产美女av久久久久小说| tocl精华| 成人永久免费在线观看视频| 国产精品久久久久久精品电影 | 伊人久久大香线蕉亚洲五| 免费看美女性在线毛片视频| 亚洲精品国产一区二区精华液| 99久久久亚洲精品蜜臀av| 黑丝袜美女国产一区| 精品久久久久久久久久免费视频| 91字幕亚洲| 亚洲国产欧美一区二区综合| 高潮久久久久久久久久久不卡| 婷婷精品国产亚洲av在线| 啪啪无遮挡十八禁网站| 看免费av毛片| 91成人精品电影| 给我免费播放毛片高清在线观看| 色播在线永久视频| 免费电影在线观看免费观看| 国产国语露脸激情在线看| 欧美国产精品va在线观看不卡| 97人妻精品一区二区三区麻豆 | 久久久久久免费高清国产稀缺| 制服诱惑二区| 极品教师在线免费播放| 不卡一级毛片| 精品第一国产精品| 黑人巨大精品欧美一区二区mp4| 国产精品二区激情视频| 欧美成人免费av一区二区三区| 国产99白浆流出| 视频区欧美日本亚洲| 久99久视频精品免费| 女性被躁到高潮视频| 国产欧美日韩精品亚洲av| 免费在线观看成人毛片| 老熟妇乱子伦视频在线观看| 天堂√8在线中文| 亚洲av电影不卡..在线观看| 动漫黄色视频在线观看| 亚洲熟妇熟女久久| 在线看三级毛片| 免费无遮挡裸体视频| 美女高潮喷水抽搐中文字幕| 最近在线观看免费完整版| 色老头精品视频在线观看| 在线天堂中文资源库| 久久久久国产一级毛片高清牌| 午夜两性在线视频| 国内毛片毛片毛片毛片毛片| 欧美另类亚洲清纯唯美| ponron亚洲| 欧美zozozo另类| 亚洲av中文字字幕乱码综合 | 99精品久久久久人妻精品| 婷婷亚洲欧美| 人人妻人人澡人人看| 两性午夜刺激爽爽歪歪视频在线观看 | 看片在线看免费视频| 国内揄拍国产精品人妻在线 | netflix在线观看网站| 亚洲精品久久成人aⅴ小说| 久99久视频精品免费| 一级a爱视频在线免费观看| 亚洲av日韩精品久久久久久密| www.精华液| 大型黄色视频在线免费观看| 91国产中文字幕| 99国产精品99久久久久| 亚洲欧美日韩无卡精品| 香蕉国产在线看| 亚洲精品一卡2卡三卡4卡5卡| 黑丝袜美女国产一区| 我的亚洲天堂| 日韩 欧美 亚洲 中文字幕| 亚洲欧美一区二区三区黑人| 中出人妻视频一区二区| 亚洲 欧美 日韩 在线 免费| 精品国产乱码久久久久久男人| 日日爽夜夜爽网站| 欧美午夜高清在线| 色综合欧美亚洲国产小说| 午夜福利在线在线| 免费看十八禁软件| 俺也久久电影网| 超碰成人久久| а√天堂www在线а√下载| 免费在线观看视频国产中文字幕亚洲| or卡值多少钱| 久久 成人 亚洲| 99在线人妻在线中文字幕| 国产aⅴ精品一区二区三区波| 欧美午夜高清在线| 级片在线观看| 国产欧美日韩一区二区三| 亚洲欧美精品综合久久99| 亚洲自拍偷在线| 99久久国产精品久久久| 悠悠久久av| 后天国语完整版免费观看| 亚洲精品在线美女| 91国产中文字幕| 亚洲成人免费电影在线观看| 日韩高清综合在线| 亚洲五月色婷婷综合| 麻豆av在线久日| 一二三四在线观看免费中文在| 一级片免费观看大全| 日韩国内少妇激情av| 国产精品精品国产色婷婷| 好看av亚洲va欧美ⅴa在| 俄罗斯特黄特色一大片| 韩国av一区二区三区四区| 日韩视频一区二区在线观看| 色播在线永久视频| 亚洲性夜色夜夜综合| 黑人操中国人逼视频| 嫩草影院精品99| 久久久久久久久中文| 天天添夜夜摸| 午夜日韩欧美国产| 一个人免费在线观看的高清视频| 麻豆一二三区av精品| www.999成人在线观看| 黄片播放在线免费| 长腿黑丝高跟| 久久青草综合色| 女生性感内裤真人,穿戴方法视频| 久久青草综合色| 亚洲国产欧洲综合997久久, | 在线观看舔阴道视频| 亚洲中文av在线| 视频在线观看一区二区三区| 亚洲第一电影网av| www日本黄色视频网| 桃色一区二区三区在线观看| 久久精品人妻少妇| 成年人黄色毛片网站| 黄色成人免费大全| 免费搜索国产男女视频| 欧美乱色亚洲激情| 制服丝袜大香蕉在线| 一区二区三区国产精品乱码| 中文字幕人妻丝袜一区二区| 色尼玛亚洲综合影院| 18禁黄网站禁片免费观看直播| 免费在线观看成人毛片| 久久欧美精品欧美久久欧美| 欧美人与性动交α欧美精品济南到| 亚洲精品一卡2卡三卡4卡5卡| 免费人成视频x8x8入口观看| 久久久久国产精品人妻aⅴ院| 嫩草影院精品99| 91老司机精品| 免费在线观看完整版高清| 国产精品久久电影中文字幕| 成人三级黄色视频| 日本撒尿小便嘘嘘汇集6| www.自偷自拍.com| 一进一出抽搐动态| 黄色 视频免费看| 久久九九热精品免费| 精品国产亚洲在线| 男人的好看免费观看在线视频 | 一本大道久久a久久精品| 国产精品爽爽va在线观看网站 | 国产99白浆流出| 日韩高清综合在线| 亚洲三区欧美一区| 1024视频免费在线观看| 亚洲精品一卡2卡三卡4卡5卡| 久久久久久久午夜电影| 看免费av毛片| 哪里可以看免费的av片| 中文字幕av电影在线播放| 变态另类丝袜制服| 99精品在免费线老司机午夜| 一二三四社区在线视频社区8| 国产又黄又爽又无遮挡在线| 亚洲成av片中文字幕在线观看| 国产片内射在线| 国产一卡二卡三卡精品| www.熟女人妻精品国产| 国产精品久久久人人做人人爽| 亚洲精品国产一区二区精华液| 日本 av在线| 亚洲精品国产一区二区精华液| 日韩 欧美 亚洲 中文字幕| 法律面前人人平等表现在哪些方面| av在线播放免费不卡| 欧美日韩亚洲国产一区二区在线观看| 免费看a级黄色片| 久久久久国产一级毛片高清牌| 男人舔女人下体高潮全视频| 日本三级黄在线观看| 两个人免费观看高清视频| 国产aⅴ精品一区二区三区波| 日本精品一区二区三区蜜桃| 久久青草综合色| 成人免费观看视频高清| 亚洲av片天天在线观看| av有码第一页| 精品卡一卡二卡四卡免费| 日韩视频一区二区在线观看| 亚洲av片天天在线观看| 久久久久久大精品| 久久久久久久久中文| 一夜夜www| 亚洲午夜精品一区,二区,三区| 夜夜爽天天搞| 黑人巨大精品欧美一区二区mp4| tocl精华| 午夜福利高清视频| 精品无人区乱码1区二区| 免费看a级黄色片| 精品日产1卡2卡| 久久精品国产亚洲av香蕉五月| 亚洲在线自拍视频| 久久久久久国产a免费观看| 久久婷婷成人综合色麻豆| 在线观看免费午夜福利视频| 亚洲av成人不卡在线观看播放网| 99热6这里只有精品| 久久人人精品亚洲av| 精品欧美一区二区三区在线| 一区福利在线观看| 伦理电影免费视频| 黄片小视频在线播放| 中文字幕最新亚洲高清| 国产精品一区二区三区四区久久 | 99国产精品一区二区三区| 啪啪无遮挡十八禁网站| 日韩中文字幕欧美一区二区| 国产av又大| 国产成人欧美| 男人的好看免费观看在线视频 | 亚洲 欧美 日韩 在线 免费| 后天国语完整版免费观看| 欧美+亚洲+日韩+国产| 国产精品亚洲av一区麻豆| 日韩欧美一区视频在线观看| 男女那种视频在线观看| 女生性感内裤真人,穿戴方法视频| 看黄色毛片网站| 51午夜福利影视在线观看| 日韩精品中文字幕看吧| 又大又爽又粗| 免费在线观看影片大全网站| 亚洲人成伊人成综合网2020| 国产成人系列免费观看| 国产激情久久老熟女| 99热只有精品国产| 脱女人内裤的视频| 国产精品精品国产色婷婷| 女人爽到高潮嗷嗷叫在线视频| av视频在线观看入口| 麻豆成人av在线观看| 国产av又大| av天堂在线播放| 一卡2卡三卡四卡精品乱码亚洲| 18禁国产床啪视频网站| 久久青草综合色| 久久中文字幕人妻熟女|