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

    Vorticity vector-potential method based on time-dependent curvilinear coordinates for two-dimensional rotating flows in closed configurations *

    2018-05-14 01:42:34YuanFuDapengZhang張大鵬XilinXie謝錫麟
    水動力學研究與進展 B輯 2018年2期
    關鍵詞:大鵬

    Yuan Fu (傅 淵), Da-peng Zhang (張大鵬), Xi-lin Xie (謝錫麟)

    Department of Aeronautics and Astronautics, Fudan University, Shanghai 200433, China

    Introduction

    The rotating driven flow is a long term concern in fluid dynamics researches which is common in nature with rich phenomenon such as geostrophic current and Ekman pumping[1,2]. The systematic studies on the topic started by Greenspan[3]are of great practical value for design and development of various kinds of fluid machinery, such as turbine machine[4-6]. In the field of spin-up problem, Wells et al.[7]studied the vortices in oscillating spin-up, and the relations between spin-up and other topics have also been investigated[8-10].

    In the studies of rotating system, the fluid governing equations are always written in the rotating frame of reference which is non-inertial. The advantage of it is the calculation domain can be fixed in the solving process. The time-dependent curvilinear coordinates also enjoy this feature[11-13]. Moreover, it is convenient to deal with the problems in irregular and varying physical configurations by using the time dependent curvilinear coordinates.

    In the field of numerical simulation, finite difference methods in the vorticity vector-potential formulation on non-staggered grids were developed[14-16]in the past decades. Vorticity-based methods enjoy advantages not only in avoiding the calculation of pressure but also in connecting with the vorticity dynamics[17]which is an important perspective on studying flow fields. It’s also convenient to deal with Coriolis term in non-inertial frame of reference[18].

    In the present study, a vorticity vector-potential method in time-dependent curvilinear coordinates is developed for two-dimensional viscous incompressible fluid in rotating system. The formula of acceleration composition in the rotating frame is derived in the curvilinear coordinates which indicates the time-dependent curvilinear coordinates is suitable for studying rotating driven flows. Tensor analysis[11,19]is adopted in deriving the component form of the governing equations. Therefore, the vorticity vector potential formulation can be discretized and solved in a regular and fixed parametric domain.

    When dealing with the boundary condition of vector-potential, the rotational vector-potential of a rigid body is adopted as the vector-potential of the fluid particles on the boundary via the non-slip condition. As to the vorticity boundary condition, the conclusion is given that in two-dimensional cases the vorticity of fluid particles at the boundary intersections is twice as the rotational angular velocity of the configuration as well as the vorticity of fluid particles at the intersection lines in three-dimensional cases. The vorticity boundary condition of the rest part can be obtained by expanding the vector-potential formula in the form of contravariant components.

    Thanks to the flexibility of time-dependent curvilinear coordinates in application, the method mentioned above is applied to triangle and curved triangle configurations respectively. The geostrophic effect, unsteady effect and curvature effect on the flow field evolutions are studied. Numerical results are presented for two-dimensional rotating driven flows.Moreover, the relationship between boundary vorticity,vorticity flux and the evolution process is preliminarily discussed.

    1. Governing equations

    1.1 Frame of reference and time-dependent curvilinear coordinates

    Time-dependent curvilinear coordinates can be introduced in a rotating system, one can set (i1,i2,i3)and (e1(t) ,e2(t) ,e3(t)) as bases of inertial frame of reference and rotating frame of reference respectively.Then the rotation of the non-inertial one can be expressed as (e1,e2,e3)(t)=(i1,i2,i3)Q(t), whereQ(t) is an orthogonal matrix. The position vector of a fluid particle in flow field can be expressed in both frames of reference asr=iX(t)=e(t)Y(t). The coordinatesY(t) in non-inertial frame can be expressed in the time-dependent curvilinear coordinates as

    wherex=x(ξ,t).xis the Eulerian coordinate, andξis the Lagrangian one.

    In the time-dependent curvilinear coordinates,the velocityVof a fluid particle in physical domain and the material derivative of a tensor fieldΦcan be respectively expressed as:

    wherer=xαgα.gαdenotes the covariant base vector of the curvilinear coordinates. Einstein summation convention is employed hereinafter.

    The velocity of a fluid particle can be decomposed into transport velocity and relative velocity as

    and the acceleration of it can be decomposed into centripetal acceleration, Coriolis acceleration and relative one

    where the relative acceleration can be expressed in the time-dependent curvilinear coordinates as

    The symboldenotes the rate of change in rotational frame of reference.

    Fig.1 Sketch of calculation of particle vorticity at vertex in a two-dimensional case

    Fig.2 Sketch of a curved triangle configuration

    Fig.3 (Color online) Sketch of the distributions in the test of the Poisson equation solver

    1.2 Vorticity vector-potential formulation

    The governing equation of vorticity in inertial frame of reference can be expressed in the time dependent curvilinear coordinates as

    where θ=V?? is the dilatation[20]. Meanwhile, in the rotating frame with the angular velocity?, the governing equation of vorticity reads

    whereωris the relative vorticity in the rotating frame of reference. If the curvilinear coordinates do not include timetexplicitly, the termwill vanish, which agrees with the governing equation in Friedlander’s book. For incompressible Newtonian fluid, dimensionless governing equation of vorticity in the inertial frame of reference reads

    where ?r/?t(x,t) is also caused by the time included in the curvilinear coordinates.

    The vorticity of a fluid particle at vertex can be determined by the rotational angular velocity of the configuration. In a two-dimensional rotating system,the vorticity of a fluid particle at the boundary intersection, namely at the vertex, is twice of the rotational angular velocity of the configuration, as shown in Fig.1.

    The calculation shows that the vorticity of a fluid particle at the vertex is only determined by the quantities of adjacent particles at the coordinate lines,namelyx1-line andx2-line, which are both on boundaries. Therefore the quantities ?V1/?x2and?V/?x1at vertex O are determined by the motion of 2 the configuration solely. The conclusion is consistent in three-dimensional cases in which the vorticity of a fluid particle at the intersection line is twice of the rotational angular velocity of the configuration.

    In the case of incompressible fluid, the vector potentialψcan be defined through the relationV=?×ψ. Then the vorticity readsω=?×(? ×ψ)= -?2ψ+ ?(? ?ψ). Under the divergence free assumption for the vector potential ??ψ=0, the governing equation of the vector-potential is ?2ψ=-ω. For two-dimensional fluid, the stream-function ψ is the degeneration of the vector-potential and the vorticity turns intoω=ωk. Here ω is a scalar function andkis the unit vector perpendicular to the configuration plane. The governing equations turn into

    2. Configuration description and boundary conditions

    2.1 Configuration description in curvilinear coordinates

    Fig.4 (Color online) Re =1000, Ω=10, t=0.05

    In a curved triangle, denoted by △ABC, the vertex A is set as the origin of the coordinates, and(ξ ,)η is set as the parametric coordinates. ξ=0 represents the vertex A, and ξ=1, η=0, η=1 represent sides BC, AB and AC respectively, as shown in Fig.2. A flat triangle configuration can be regarded as a special case of it.

    The position vector of a point P in the curved triangle configuration can be expressed as

    whererb(η) is a vector from the vertex A to side BC through the point P. The vector valued mappingr( ξ ,η) establishes a diffeomorphism between the parametric domainDξηand the interior of the configuration ΔABC. The covariant basis vectors of the coordinates (ξ ,)η can be calculated as

    and the covariant components of the metric tensor are

    The determinant of metric matrix is

    Fig.5 (Color online) Re =1000, Ω=10, t=0.79

    The contravariant basis vectors of the coordinates(ξ ,)η are

    The non-zero second kind of Christoffel symbols can be calculated as:

    and the non-zero first kind of Christoffel symbols can be calculated as:

    In the case of flat triangle configurations, we have

    Fig.6 (Color online) Re =1000, Ω=10, t=2.36

    While in the case of curved triangle configurations,the curve of BC boundary is set as

    which is determined by two parameters λ and μ.Here λ determines the amplitude of the bend and μ determines the wave number of it. sinμ is forced to be zero, namely μ=kπ, in order to keep the end of vectorrCto coincide with the vertex C of BC boundary.

    2.2 Boundary conditions of vector-potential and vorticity

    Fig.7 Sketch of the velocity distribution of clockwise circulation in the flow field

    The velocity of a fluid particle on the boundary of the configuration equals to the velocity of the rigid boundary point. So the rotational vector-potential of the rigid body is adopted as the vector-potential of fluid particles on the boundary, namely

    Fig.8 (Color online) Re =1000, Ω=10, t=19.01

    In the case of two-dimensional incompressible flow,the relation between stream-function ψ and contravariant components of velocity can be expressed in the parametric coordinates (ξ ,)η as:

    wheregis the determinant of the coordinates.

    Due to the relation2ψ=ω-? , the boundary condition of vorticity can be calculated by the expansion of the stream-function

    In the case of three-dimensional flow, the derivatives of the contravariant components of the vector potential can be generally expressed as:

    In the case of two-dimensional flow, it turns to be

    It’s also worth pointing out that the singularity at the vertex A included in the formula mentioned above can be overcome by using the result that the vorticity on the vertex is twice of the rotational angular velocity of the configuration as mentioned in Section 1.2.

    3. Numerical method and validation

    3.1 Numerical scheme

    A finite difference method is employed to carry out the spatial discretization. To be specific, a nonisometric-interval Lagrangian five-point interpolation method is adopted. The calculation grid number is 80×60. Meanwhile, the prediction-correction method is adopted to discretize time, and the time step is set as dt=1× 10-5.

    Table 1 Maximum errors of geometry quantities in the curved triangle configuration with parameters λ=0.01, μ=5π

    In the test of grid quality, the geometry quantities including components of metric tensor and two kinds of Christoffel symbols are numerically calculated. The results are compared with analytical solutions. In the case of curved triangle configuration with parameters λ= 0.01, μ=5π, the numerical errors are reported in the Table 1. The symbol “NaN” implies a zero solution of the geometry quantity. The data in the table shows that the magnitude of relative errors on the whole domain is no larger than 10-4, except the data 3.6% caused by the almost-zero analytical solutions. It indicates that the grid mentioned above is acceptable.

    Fig.9 (Color online) Re =1000, Ω=1, t=100.22

    3.2 Validation of the Poisson equation solver

    The governing equation of the stream-function is a scalar Poisson equation ?2ψ / ?x2+ ?2ψ / ?y2+f(x,y) =0. To test the solver, the source term is set asf(x,y)=-4 . The function ψ(x,y)=x2+y2is supposed to be the exact solution of the Poisson equation and the boundary condition is given correspondingly.Then the Poisson equation is mapped to the parametric domain to iterate for a numerical solution. The distributions of the numerical solution and absolute error in the curved triangle configuration are shown in Fig.3.

    Fig.10 (Color online) Re =2000, Ω=10, t=10.05

    The result shows that the maximum value of absolute error on the whole domain is 7.58×10-5, and the maximum value of relative error is 4.70×10-4. It also indicates that the numerical method mentioned above is suitable for accurately simulating the rotating system.

    4. Numerical results and discussion

    4.1 Constant angular velocity cases with flat boundary

    In cases where angular velocity is constant, the flow field evolution can be divided into several stages.Take the case of angular velocity Ω=10,Re=1000 as an example, where the characteristic lengthLinRe=ρUL/μ is set as the length of the triangle side. At the beginning, high vorticity is generated in the middle of the boundaries by a sudden start of the configuration rotation, as shown in Fig.4.

    Fig.11 (Color online) Re =1000, Ω=5.5+4.5sin(5.5t),t=0.45

    In the graphs here and after, thex-axis measures the scaling length of the boundary, they- axis on the left and right measure the corresponding vorticity flux and vorticity respectively, where the vorticity flux ?/ω?? on the boundary represents the intensity of vorticity transportation. ? is the coordinate in direction ofnwhich is perpendicular to the boundary in configuration plane. When the boundary vorticity has the same sign as vorticity flux, the vorticity is enhancing, and vice versa.

    Fig.12 (Color online) Re =1000, Ω=5.5+4.5sin(5.5t),t=0.99

    At the sudden start, the vorticity at vertex B and C isas mentioned in Section1.2.The vorticity at the boundary AB and AC is positive, as shown in Fig.4. According to the Rolle theorem,the vorticity has a maximum value at the BC boundary. Actually, it reaches the peak value at the middle point of BC due to the symmetry. Considering the stress tensor at the boundary,t=-pI+2μD, whereDis the rate of strain tensor, the stresst?n=-pn+μω×nacts on BC. The shear stress μωis a resistance against configuration rotation. The reactive force of it is against the clockwise circulation, which is dominated by the geostrophic effect. In this case the geostrophic effect is stronger than the effect of shear stress.

    Att=0.79, one end of the region where vorticity is concentrated departs from the boundary to the interior of the flow field. A clockwise circulation is formed in the middle of it, while three small counter-clockwise circulations are formed in the corner regions, as shown in Fig.5. There is an opposite tendency between the distributions of vorticity and vorticity flux at the boundary AB and AC. At boundary BC, there is an about π/2 phase difference between them.

    Fig.13 Re =1000, Ω=5.5+4.5sin(5.5t), vorticity and vorticity flux at boundary AB

    Then, the regions of vorticity concentration are connected with each other to form a loop. The flow field is dominated by the single clockwise circulation,as shown in Fig.6. Compared with the distribution of boundary vorticity at the momentt=0.79, the vorticity flux at boundary reduces it by transporting it to inside the flow field. The amplitude of both boundary vorticity and vorticity flux decreases.

    The dimensionless momentum equation in the rotating frame of reference reads

    whereRo=U/ΩLandEc= ν/ΩL2denote Rossby number and Ekman number respectively.When the Rossby number is small, the convective term can be omitted. The following equation indicates the geostrophic equilibrium

    The quasi-geostrophic clockwise flow in the middle of the flow field indicates that the ?ppoints to the center of it, as shown in Fig.7, whererCis a vector from the vertex A to the triangle center, anderis the normalization of the position vectorrwhich also starts from the vertex A.

    Fig.14 (Color online) Re =1000, Ω=5.5+4.5sin(5.5t),t =1.20,

    Whenthere is. The result isp(rC)>p(r). On the other hand, when, there isp(r) -p(rC). The result is alsop(rC)>p(r). In conclusion,p(rC) is the maximum value ofpin the flow field. The direction of the circulation agrees with the result of geostrophic balance.

    Att=19.01, the loop of vorticity concentration begins to dissipate.The value of vorticity flux at the boundaries is near to zero. However, the non-zero vorticity is maintained at the boundaries, as shown in Fig.8. In the process of the flow field evolution, the vorticity flux always suppresses the boundary vorticity. Then the boundary vorticity is reduced continually while the vorticity flux is near to zero on AB and AC boundaries. The vorticity flux at BC boundary is maintained a symmetric distribution and the state of flow field approaches to a geostrophic equilibrium.

    Fig.15 (Color online) Re =1000, Ω=5.5+4.5sin(5.5t),t=1.62

    In the cases of constant rotational angular velocity of configuration, such as (Re=1000, Ω=1),(Re=2000, Ω=1), and (Re=2000, Ω=10), the flow field evolution goes through similar stages.. The dominated circulation in flow field is always in the opposite direction of the configuration rotation, and the amplitude of boundary vorticity and vorticity flux decreases in the evolution process. Especially, in case of (Re=1000, Ω=1), the final state of flow field is like the rotation of a rigid body. The concentration of vorticity vanishes, and the amplitude of the relative velocity drops to zero almost everywhere in the flow field, as shown in Fig.9.

    Fig.16 (Color online) μ=3π, λ=0.05, Re =2000,Ω =10, t=0.47

    With a larger Reynolds number, there appears a centrifugal effect in the flow field, so that the flow pattern is slightly twisted. As shown in Fig.10, in the case of (Re=2000, Ω=10), the region of high vorticity value concentration in the middle of the flow field has an elliptical outline rather than a circular one.

    4.2 Varying angular velocity cases with flat boundary

    In a rotating system, rich phenomena can be created by the unsteady effect. In the case ofRe=1000,Ω =5.5+4.5sin(5.5t), the evolution process of the flow field is described as follows.

    Fig.17 (Color online) μ=3π, λ=0.05, Re =2000,Ω =10, t=17.28

    Firstly, as in the case of constant rotational angular velocity of configuration, the high vorticity value is generated in the middle of the boundaries and enters the flow field. A clockwise circulation is formed in the middle of the flow field and three small counter-clockwise circulations are formed at the corners, as shown in Fig.11. Due to the unsteady effect, with the decrease of the rotational angular velocity of the configuration, the counter-clockwise circulations are enhanced to dominant the whole flow field rather than vanish in the evolution. Then the counter-clockwise circulations touch each other and begin to merge.

    Att=0.99, the counter-clockwise circulations merge in the middle of the flow field. A unified counter-clockwise circulation is formed, as shown in Fig.12, which is not found in the constant rotational angular velocity cases. The unsteady effect could change the flow field significantly.

    Fig.18 μ=3π , λ=0.05, Re =2000, Ω=10, vorticity and vorticity flux at boundary BC and AC

    At boundary AB, there is an opposite tendency between the distributions of boundary vorticity and vorticity flux. Fromt=0.63 tot=0.99, the position of the peak value moves from 0.2 to 0.5 with the motion of the counter-clockwise circulations, as shown in Fig.13. In this process, the amplitudes of both the boundary vorticity and vorticity flux decrease.The counter-clockwise circulations move from the corner regions to the middle of the flow field.

    With the increase of the rotational angular velocity of the configuration, the new clockwise circulations are formed and enhanced in the corner regions. Att=1.20, there is a counter-clockwise circulation in the interior region of the flow field and a clockwise one in the outer region of it, as shown in Fig.14.

    Fig.19 (Color on line) μ=7π, λ= 0.01, Re =2000,Ω =10,t=11.31

    The acceleration of a fluid particle is

    where the term

    represents the unsteady effect. In the case of constant angular velocity rotation, the direction of the dominated circulation in the flow field is always clockwise which is opposite to the configuration rotation, and the geostrophic balance equation is expressed as Eq. (20). In the varying rotational angular velocity cases, the equation turns to

    Fig.20 (Color online) μ=7π, λ=0.01, Re =2000,Ω =10, t=24.50

    With the increase of the rotational angular velocity of the configuration, namelythe clockwise circulation is enhanced, and with the decrease of it, namelythe counter-clockwise one is enhanced instead.

    With the periodic change of the rotational angular velocity of the configuration, the multi-vortex structure appears in the flow field as shown in Fig.15.Att=1.62, the flow field is dominated by three clockwise circulations, and there are small counterclockwise circulations near the boundary at the same time. The flow field evolution is quasi-periodic. Att=2 .55 there is a clock wise circulation in the flow field which is similar to the initial state of the flow field.

    Fig.21 (Color online) μ=3π, λ=0.05, Re =2000,Ω=5.5+4.5sin(5.5t), t=0.39

    In the cases of varying rotational angular velocity with flat boundary, such as (Re=1000,Ω=15+5sin(15t)), (Re=2000,Ω =5.5+4.5sin(5.5t)), and(Re=2000, Ω=15+5sin(15t)), the flow field evolution goes through similar stages as discussed above.

    4.3 Constant angular velocity cases with curved boundary

    Besides the unsteady effect, the curvature effect on the flow field evolution can not be omitted. In the case ofRe=2000, Ω=10 with BC boundary parameters μ=3π, λ=0.05, the evolution process of the flow field is described as follows.

    In the initial stage, the dominated clockwise circulation is formed in the flow field. Meanwhile there is a small counter-clockwise circulation on the BC boundary as in backward-facing step flow. The counter-clockwise circulation enhances the dominated one. Moreover, the dominated clockwise circulation enhances another counter-clockwise circulation in the region of corner A, as shown in Fig.16.

    Then the small counter-clockwise circulations slowly vanish, and the dominated clockwise circulation is fixed on the bend of BC, rotating with the configuration, as shown in Fig.17. Compared with the dominated circulations in the flat boundary configurations, the dominated one in this case shares the same direction with them which is determined by the geostrophic equilibrium. The curvature effect merely changes the position of it in this case.

    There is an about π/2 phase difference between vorticity and vorticity flux on boundary BC. Att=0.47, there is an opposite tendency between the distributions of vorticity and vorticity flux at the boundary AC. The peak values of them appear at about 0.4 position, where the clockwise circulation touches the boundary. Att=17.28, the peak value position moves to 0.5, both vorticity and vorticity flux are decreased, as shown in Fig.18. In the case ofRe=2000, Ω=10 with boundary parameters μ=5π, λ=0.02, the flow field evolution goes through a similar process. The dominated clockwise circulation is fixed on the bend of BC eventually.

    However, in the case of curved boundary with parameters μ=7π, λ=0.01, the flow field evolution is changed significantly by the curvature effect.After the formation of clockwise circulation in the middle of the flow field, the counter-clockwise circulation at region of corner C is enhanced rather than suppressed.

    Att=11.31, there are two counter-rotating circulations in the flow field. They are both at the bends of boundary BC, as shown in Fig.19, which is not seen in other cases with the same Reynolds number and rotational angular velocity of configuration.

    Eventually the counter-clockwise circulation slowly vanishes, as shown in Fig.20. The final state of the flow field is still a geostrophic equilibrium.However, the size of the clockwise circulation in this case is smaller than in others as the result of curvature effect.

    4.4 Varying angular velocity cases with curved boundary

    Unsteady effect and curvature effect work together in the following case. The parameters of it areRe=2000, Ω=5.5+4.5sin(5.5t), and the parameters of BC boundary are μ=3π, λ=0.05.

    At the initial phase, the dominated clockwise circulation is formed in the middle of the flow field.Two small counter-clockwise circulations are formed on the BC boundary and in the region of corner A respectively, as shown in Fig.21. With the decrease of the rotational angular velocity of the configuration,the counter-clockwise circulations are enhanced due to the unsteady effect.

    Fig.22 (Color online) μ=3π, λ=0.05, Re =2000,Ω=5.5+4.5sin(5.5t), t=0.78

    In the process of the counter-clockwise circulations enhancing, they are gathered close to each other along the boundary AB, as shown in Fig.22. The curvature effect breaks the rotational symmetry of the motion path. Att=0.99, the counter-clockwise circulations merge in the middle of the flow field. A unified counter-clockwise circulation with new clockwise circulations are formed. With the increase of the rotational angular velocity of the configuration,the clockwise circulations are enhanced. Fromt=0.99 tot=1.35, the rotational angular velocity of the configuration keeps increasing, the clockwise circulations are enhanced due to the unsteady effect.Eventually a unified clockwise circulation is formed in the flow field.

    The vorticity flux at the curved boundary BC can be decomposed into

    here ? is the coordinate in direction ofn,Vη is the component of velocity in direction ofτ, which is the unit tangent vector of BC, and ληis the curvature of it. The distributions of these terms at different moments are shown in Fig.23. The term-ληω contains the curvature effect at the boundary BC, which vanishes in case of flat boundary configuration. Around every peak value of vorticity flux,curvature term plays an important role in the composition of it. There is an about π/2 phase difference between the boundary vorticity and vorticity flux. Fromt=0.39 tot=0.78, with the decrease of the rotational angular velocity of the configuration, the amplitudes of both vorticity and vorticity flux on BC decrease, and vice versa fromt=0.99 tot=1.35.

    Fig.23 μ=3π,λ= 0.05, Re =2000, Ω =5.5+4.5sin(5.5t), vorticity and vorticity flux at boundary BC

    In this periodic process, the increase and decrease of the rotational angular velocity of the configuration cause the enhancing of clockwise circulations and counter-clockwise ones respectively. The flow field evolution is quasi-periodic.

    5. Conclusions

    In this study, a numerical method based on vorticity vector-potential formulation is developed in the time-dependent curvilinear coordinates. The composition formulas of velocity and acceleration in rotating frame of reference are included in the curvilinear coordinates description. The centripetal acceleration, Coriolis acceleration and relative one can be simply expressed as two terms in this coordinates.The parametric domain is fixed and regular in the calculation process. The method is applied to the simulations of two-dimensional viscous incompressible rotating driven flows in closed configurations by overcoming the singularity at rotation axis.The numerical results in triangle and curved triangle configurations are presented by demonstrating the distributions of relative stream-function, relative vorticity and relative velocity. The flow field evolutions are discussed and classified by the shape of the configurations as well as the rotational angular velocity of them.

    In the case of constant rotational angular velocity,high value of vorticity appears at the boundaries due to the sudden start of rotation. The dominated circulation which is determined by the geostrophic effect in the flow field is always in opposite direction of the configuration rotation. At AB and AC boundaries, there is an opposite tendency between the distributions of vorticity and vorticity flux, and an about π/2 phase difference between them at BC boundary. The flow fields approach geostrophic equilibrium in the final stage of the evolution. In the case of low angular velocity with small Reynolds number,the flow field behaves like a rotating rigid body eventually. Curvature effect plays different roles in different cases. Usually it changes flow fields on details, such as changing the position of circulations and the route of their motions. However, it can also change the evolution significantly, such as in the case of (μ=7π ,λ=0.01,Re=2000,Ω=10). The curvature effect doesn’t influence the direction of the dominated circulation in a flow field.

    In the case of varying rotational angular velocity of the configuration, the evolution is quasi-periodic because of the periodic change of rotational angular velocity. When the rotational angular velocity of the configuration increases, the clockwise circulations in the flow field are enhanced and vice versa. The unsteady effect changes the direction of the dominated circulation from time to time, and the flow field is much more complicated than in case of constant rotational angular velocity. However, no matter how complex the evolution is, the pattern of it varies cyclically, and the pattern similar to the initial state can always be found.

    The proposed method can be extended to three-dimensional situations with even more complicated configurations. It will be useful to calculate flow fields in different containers such as engines or rotating pumps. The extensions will be considered in further researches.

    [1] Hardenberg J., McWilliams J. C., Provenzale A. et al.Vortex merging in quasi-geostrophic [J].Journal of Fluid Mechanics, 2000, 412: 331-353.

    [2] Plumley M., Julien K., Marti P. et al. The effects of Ekman pumping on quasi-geostrophic Rayleigh-Benard convection [J].Journal of Fluid Mechanics, 2016, 803:51-71.

    [3] Greenspan H. P. The theory of rotating fluids [M].Cambridge, UK: Cambridge University Press, 1968.

    [4] You D., Wang M., Moin P. et al. Large-eddy simulation analysis of mechanisms for viscous losses in a turbomachinery tip-clearance flow [J].Journal of Fluid Mechanics,2007, 586: 177-204.

    [5] Li L. J., Zhou S. J. Numerical simulation of hydrodynamic performance of blade position-variable hydraulic turbine[J].Journal of Hydrodynamics, 2017, 29(2): 314-321.

    [6] Zhou P. J., Wang F. J., Yang Z. J. et al. Investigation of rotating stall for a centrifugal pump impeller using various SGS models [J].Journal of Hydrodynamics, 2017, 29(2):235-242.

    [7] Wells M. G., Clercx H. J. H., Van Heijst G. J. F. Vortices in oscillating spin-up [J].Journal of Fluid Mechanics,2007, 573: 339-369.

    [8] Foster M. R., Munro R. J. The linear spin-up of a stratified,rotating fluid in a square cylinder [J].Journal of Fluid Mechanics, 2012, 712: 7-40.

    [9] Flor J. B., Ungarish M., Bush J. W. M. Spin-up from rest in a stratified fluid: boundary flows [J].Journal of Fluid Mechanics, 2002, 472: 51-82.

    [10] Suh Y. K., Choi Y. H. Study on the spin-up of fluid in a rectangular container using Ekman pumping models [J].Journal of Fluid Mechanics, 2002, 458: 103-132.

    [11] Xie X. L. Modern Tensor analysis with applications in continuum mechanics [M]. Shanghai, China: Fudan University Press, 2014(in Chinese).

    [12] Xie X. L. On two kinds of differential operators on general smooth surfaces [J].Journal of Fudan University (Natural Science), 2013, 52(5): 688-711.

    [13] Xie X. L., Chen Y., Shi Q. Some studies on mechanics of continuous mediums viewed as differential manifolds [J].Science China Physics, Mechanics and Astronomy, 2013,56(2): 432-456.

    [14] E W., Liu J. Finite difference methods for 3D viscous incompressible flows in the vorticity vector potential formulation on non-staggered grids [J].Journal of Computational Physics, 1997, 138(1): 57-82.

    [15] Meitz H. L., Fasel H. F. A compact-difference scheme for the Navier-Stokes equations in vorticity-velocity formulation [J].Journal of Computational Physics, 2000, 157(1):371-403.

    [16] Olshanskii M. A., Rebholz L. G. Velocity-vorticity-heicity formulation and a solver for the Navier-Stokes equations[J].Journal of Computational Physics, 2010, 229(11):4291-4303.

    [17] Wu J. Z., Ma H. Y., Zhou M. D. Vorticity and vortex dynamics [M]. Berlin, Heidelberg, Germany: Springer,2005.

    [18] Luo H., Bewley T. R. On the contravariant form of the Navier-Stokes equations in time-dependent curvilinear coordinates systems [J].Journal of Computational Physics,2004, 199(1): 355-375.

    [19] Dubrovin B. A., Fomenko A. T., Novikov S. P. Modern geometry-methods and applications [M]. Berlin,Heidelberg, Germany: Springr-Verlag, 2011.

    [20] Chen Y., Xie X. Vorticity vector-potential method for 3D viscous incompressible flows in time-dependent curvilinear coordinates [J].Journal of Computational Physics,2016, 312: 50-81.

    猜你喜歡
    大鵬
    周鵬飛:大鵬展翅 跨界高飛
    華人時刊(2022年7期)2022-06-05 07:33:46
    Beating standard quantum limit via two-axis magnetic susceptibility measurement
    看圖紙
    當代工人(2019年20期)2019-12-13 08:26:11
    三棱錐中的一個不等式
    Proton Beam Generated by Multi-Lasers Interaction with Rear-Holed Target
    李大鵬:打造縱向、橫向全域發(fā)展的蘇交科
    中國公路(2017年14期)2017-09-26 11:51:42
    《縫紉機樂隊》懟大鵬發(fā)布會 娜扎奪大鵬銀幕初吻
    劉業(yè)偉、王大鵬設計作品
    藝術評論(2017年3期)2017-05-04 12:49:00
    非誠勿擾
    AComparingandContrastingAnalysisofCooperationandPoliteness
    欧美中文综合在线视频| 日韩精品青青久久久久久| 岛国在线观看网站| 欧美绝顶高潮抽搐喷水| 激情在线观看视频在线高清| 久久久精品大字幕| 男人和女人高潮做爰伦理| 97超级碰碰碰精品色视频在线观看| 日韩欧美在线二视频| 欧美最黄视频在线播放免费| 国内精品久久久久精免费| av专区在线播放| 高清在线国产一区| 真实男女啪啪啪动态图| 国产久久久一区二区三区| 欧洲精品卡2卡3卡4卡5卡区| 欧洲精品卡2卡3卡4卡5卡区| 亚洲天堂国产精品一区在线| 亚洲国产欧美网| 内射极品少妇av片p| 日本熟妇午夜| 亚洲最大成人中文| 夜夜看夜夜爽夜夜摸| 97碰自拍视频| 亚洲美女黄片视频| 丁香欧美五月| 3wmmmm亚洲av在线观看| 99久久精品热视频| 午夜福利免费观看在线| 亚洲精华国产精华精| 18禁黄网站禁片午夜丰满| xxxwww97欧美| 色尼玛亚洲综合影院| 天堂网av新在线| netflix在线观看网站| 一级毛片女人18水好多| 国产探花在线观看一区二区| 免费一级毛片在线播放高清视频| 国产探花在线观看一区二区| 看黄色毛片网站| 中文资源天堂在线| 欧美日韩瑟瑟在线播放| 国产乱人伦免费视频| 白带黄色成豆腐渣| 国产午夜福利久久久久久| 日本五十路高清| 国产一区二区亚洲精品在线观看| 最新在线观看一区二区三区| 日韩精品中文字幕看吧| 精华霜和精华液先用哪个| 一个人免费在线观看电影| av片东京热男人的天堂| 757午夜福利合集在线观看| 村上凉子中文字幕在线| 欧美日韩乱码在线| 黄色丝袜av网址大全| 亚洲专区国产一区二区| 中文字幕人妻丝袜一区二区| 一本一本综合久久| 国产精品电影一区二区三区| 亚洲精品一区av在线观看| 久久久精品欧美日韩精品| 热99re8久久精品国产| 99久久精品一区二区三区| 12—13女人毛片做爰片一| 老汉色av国产亚洲站长工具| 午夜激情福利司机影院| 欧美+日韩+精品| 午夜福利18| 一个人免费在线观看电影| 蜜桃亚洲精品一区二区三区| av福利片在线观看| 欧美黑人欧美精品刺激| 亚洲国产色片| 欧美乱码精品一区二区三区| 午夜影院日韩av| 国内毛片毛片毛片毛片毛片| 亚洲国产精品久久男人天堂| 12—13女人毛片做爰片一| 国产三级黄色录像| 欧美成人免费av一区二区三区| 国内精品一区二区在线观看| 黄色女人牲交| 精品一区二区三区视频在线 | 久久久久久九九精品二区国产| 亚洲精品粉嫩美女一区| 国产亚洲欧美98| 两性午夜刺激爽爽歪歪视频在线观看| 婷婷六月久久综合丁香| 精品一区二区三区视频在线 | 一夜夜www| 精品日产1卡2卡| h日本视频在线播放| 久久婷婷人人爽人人干人人爱| 亚洲精品影视一区二区三区av| 国产久久久一区二区三区| 九九在线视频观看精品| 色综合站精品国产| 欧美一级毛片孕妇| 女同久久另类99精品国产91| 国产欧美日韩精品亚洲av| 麻豆成人午夜福利视频| 久久久久免费精品人妻一区二区| 国产精品久久久久久精品电影| 国产精品 国内视频| 最新在线观看一区二区三区| 国产欧美日韩精品一区二区| 成年版毛片免费区| 精品电影一区二区在线| 久久久久久久久大av| 啪啪无遮挡十八禁网站| 欧美色欧美亚洲另类二区| 最近最新中文字幕大全电影3| 国产熟女xx| 日本熟妇午夜| 亚洲成人久久性| 亚洲美女视频黄频| 国产一区二区激情短视频| 国产黄色小视频在线观看| 亚洲人成伊人成综合网2020| 一级毛片女人18水好多| 高潮久久久久久久久久久不卡| 精品国内亚洲2022精品成人| 丁香六月欧美| 欧美绝顶高潮抽搐喷水| netflix在线观看网站| 久久香蕉精品热| 午夜影院日韩av| 男女床上黄色一级片免费看| 香蕉丝袜av| 国产亚洲精品久久久com| 久久久国产成人精品二区| 欧美一级a爱片免费观看看| 岛国在线观看网站| 国内少妇人妻偷人精品xxx网站| 国产欧美日韩一区二区三| 欧美成人a在线观看| 欧美一级毛片孕妇| 久久精品国产综合久久久| 91久久精品电影网| 精品国产美女av久久久久小说| 香蕉丝袜av| 亚洲乱码一区二区免费版| 高清日韩中文字幕在线| 99riav亚洲国产免费| 日日干狠狠操夜夜爽| 老司机福利观看| 一个人看视频在线观看www免费 | 久久6这里有精品| 女警被强在线播放| 精品人妻一区二区三区麻豆 | 免费av观看视频| 精品国产亚洲在线| 国产不卡一卡二| 黄色视频,在线免费观看| 精品熟女少妇八av免费久了| 久久人妻av系列| 欧美日韩综合久久久久久 | 婷婷丁香在线五月| 国产伦在线观看视频一区| 亚洲精品亚洲一区二区| 有码 亚洲区| 一本久久中文字幕| 国产午夜精品论理片| 国产高清激情床上av| 一本一本综合久久| 精品福利观看| 淫秽高清视频在线观看| 亚洲天堂国产精品一区在线| 最好的美女福利视频网| 久久性视频一级片| 国产精品亚洲美女久久久| 婷婷六月久久综合丁香| 午夜激情福利司机影院| 亚洲国产高清在线一区二区三| 琪琪午夜伦伦电影理论片6080| 亚洲第一欧美日韩一区二区三区| 老熟妇仑乱视频hdxx| 国产高清有码在线观看视频| 天堂av国产一区二区熟女人妻| eeuss影院久久| 亚洲最大成人手机在线| 国产伦在线观看视频一区| 香蕉av资源在线| 欧美日韩福利视频一区二区| 一个人看的www免费观看视频| 亚洲精品美女久久久久99蜜臀| 欧美丝袜亚洲另类 | 国产伦在线观看视频一区| 国产精品久久久人人做人人爽| 亚洲人与动物交配视频| 在线观看免费视频日本深夜| a级毛片a级免费在线| 欧美极品一区二区三区四区| 精品一区二区三区人妻视频| 母亲3免费完整高清在线观看| 窝窝影院91人妻| 97人妻精品一区二区三区麻豆| 日韩高清综合在线| 成人高潮视频无遮挡免费网站| 亚洲熟妇中文字幕五十中出| 亚洲国产精品久久男人天堂| 无遮挡黄片免费观看| 十八禁人妻一区二区| 9191精品国产免费久久| 色综合亚洲欧美另类图片| 少妇人妻精品综合一区二区 | 亚洲 国产 在线| 欧美黑人巨大hd| 精品久久久久久久人妻蜜臀av| 一a级毛片在线观看| 最近在线观看免费完整版| www日本在线高清视频| 男女床上黄色一级片免费看| 成人18禁在线播放| 少妇人妻一区二区三区视频| 国产麻豆成人av免费视频| 手机成人av网站| 亚洲熟妇熟女久久| 校园春色视频在线观看| 午夜影院日韩av| 亚洲最大成人手机在线| 欧美最新免费一区二区三区 | 嫩草影院入口| 香蕉久久夜色| 97超视频在线观看视频| 最新在线观看一区二区三区| 两个人的视频大全免费| 中亚洲国语对白在线视频| 欧美高清成人免费视频www| 九九久久精品国产亚洲av麻豆| 18禁黄网站禁片午夜丰满| 亚洲专区国产一区二区| 最近最新中文字幕大全免费视频| 中文资源天堂在线| 久久久久久国产a免费观看| 国产成人啪精品午夜网站| av在线蜜桃| 亚洲欧美日韩卡通动漫| 久久精品国产综合久久久| 97超视频在线观看视频| 高清在线国产一区| 午夜a级毛片| 中文字幕精品亚洲无线码一区| 变态另类成人亚洲欧美熟女| 亚洲七黄色美女视频| 美女被艹到高潮喷水动态| 黄色女人牲交| 亚洲国产中文字幕在线视频| 精品久久久久久,| av女优亚洲男人天堂| 国产精品一区二区三区四区免费观看 | 免费看美女性在线毛片视频| 精品国内亚洲2022精品成人| 亚洲av免费在线观看| 婷婷丁香在线五月| 亚洲欧美日韩无卡精品| 99久久精品国产亚洲精品| 色哟哟哟哟哟哟| 免费人成视频x8x8入口观看| 久久精品国产自在天天线| 国产成人av教育| 在线观看免费视频日本深夜| 亚洲精品影视一区二区三区av| 午夜日韩欧美国产| 3wmmmm亚洲av在线观看| 十八禁人妻一区二区| 怎么达到女性高潮| 成人av一区二区三区在线看| 亚洲18禁久久av| 桃红色精品国产亚洲av| 亚洲天堂国产精品一区在线| 人妻夜夜爽99麻豆av| 露出奶头的视频| www国产在线视频色| 黄色成人免费大全| 中文字幕人妻丝袜一区二区| 日本撒尿小便嘘嘘汇集6| 免费看光身美女| 18禁国产床啪视频网站| 麻豆久久精品国产亚洲av| 搡女人真爽免费视频火全软件 | 99久久精品热视频| 最后的刺客免费高清国语| 成人精品一区二区免费| 2021天堂中文幕一二区在线观| 性色avwww在线观看| 久久精品国产综合久久久| 欧美午夜高清在线| 亚洲七黄色美女视频| 国产高清三级在线| 亚洲真实伦在线观看| 国产高清有码在线观看视频| 一级毛片女人18水好多| 悠悠久久av| 香蕉丝袜av| 丰满人妻熟妇乱又伦精品不卡| 欧美区成人在线视频| 中文字幕人妻丝袜一区二区| 欧美丝袜亚洲另类 | 国产精品嫩草影院av在线观看 | 中文字幕人妻丝袜一区二区| 国内精品久久久久久久电影| 欧美绝顶高潮抽搐喷水| 日韩 欧美 亚洲 中文字幕| a在线观看视频网站| 精品国内亚洲2022精品成人| 婷婷精品国产亚洲av在线| 亚洲av熟女| 国产精品一区二区免费欧美| 久久久久久九九精品二区国产| 日本免费a在线| 18禁黄网站禁片午夜丰满| 老鸭窝网址在线观看| 麻豆一二三区av精品| 亚洲精品成人久久久久久| 欧美成人a在线观看| 日日夜夜操网爽| 日韩成人在线观看一区二区三区| 日本在线视频免费播放| 在线播放无遮挡| 国产免费一级a男人的天堂| 亚洲欧美一区二区三区黑人| 夜夜躁狠狠躁天天躁| 少妇丰满av| 熟女人妻精品中文字幕| 婷婷精品国产亚洲av| 久久这里只有精品中国| 日韩高清综合在线| 天天一区二区日本电影三级| 十八禁网站免费在线| 日本撒尿小便嘘嘘汇集6| 白带黄色成豆腐渣| 国产亚洲欧美在线一区二区| 久久草成人影院| 一区二区三区国产精品乱码| 国产激情偷乱视频一区二区| 成人国产一区最新在线观看| 一个人观看的视频www高清免费观看| 精品国内亚洲2022精品成人| 欧美一区二区精品小视频在线| 丰满人妻一区二区三区视频av | 一本一本综合久久| 搡老熟女国产l中国老女人| 国产精品一及| 免费在线观看成人毛片| 欧美不卡视频在线免费观看| 国产野战对白在线观看| 天美传媒精品一区二区| 成年免费大片在线观看| 国产熟女xx| 久久精品国产亚洲av香蕉五月| 欧美日韩乱码在线| 日韩大尺度精品在线看网址| 美女免费视频网站| 18禁裸乳无遮挡免费网站照片| 日本三级黄在线观看| 亚洲色图av天堂| 丰满的人妻完整版| 国产真实乱freesex| 18+在线观看网站| 午夜福利欧美成人| 精品熟女少妇八av免费久了| 白带黄色成豆腐渣| 九九热线精品视视频播放| 国产成人啪精品午夜网站| 日韩免费av在线播放| 一边摸一边抽搐一进一小说| 亚洲国产日韩欧美精品在线观看 | 亚洲国产中文字幕在线视频| 国产爱豆传媒在线观看| 欧美zozozo另类| 啦啦啦韩国在线观看视频| 欧美日韩福利视频一区二区| 久久久久久久久中文| 一进一出好大好爽视频| 在线看三级毛片| 久久久久久人人人人人| 一本久久中文字幕| 99久久九九国产精品国产免费| 露出奶头的视频| 国产色爽女视频免费观看| 日本五十路高清| 一区二区三区免费毛片| 波多野结衣巨乳人妻| 人人妻,人人澡人人爽秒播| 欧美日韩瑟瑟在线播放| 黄片小视频在线播放| 午夜精品在线福利| 国产精品一区二区三区四区免费观看 | 免费在线观看影片大全网站| 最近最新中文字幕大全免费视频| 国产激情偷乱视频一区二区| 亚洲美女视频黄频| aaaaa片日本免费| 国产极品精品免费视频能看的| 亚洲av熟女| 免费av观看视频| 国产国拍精品亚洲av在线观看 | 精品国产超薄肉色丝袜足j| 国产伦人伦偷精品视频| 日日夜夜操网爽| 久久久久久久久大av| 特大巨黑吊av在线直播| 国产乱人视频| 十八禁人妻一区二区| 婷婷亚洲欧美| 国产精品女同一区二区软件 | www.www免费av| 精品福利观看| 亚洲最大成人手机在线| 国产三级中文精品| 又爽又黄无遮挡网站| 色噜噜av男人的天堂激情| 国产精品永久免费网站| 五月玫瑰六月丁香| 午夜视频国产福利| 亚洲av第一区精品v没综合| 免费无遮挡裸体视频| 欧美成人一区二区免费高清观看| 国产极品精品免费视频能看的| 成人av在线播放网站| 男人舔奶头视频| 又爽又黄无遮挡网站| 午夜a级毛片| 亚洲天堂国产精品一区在线| 亚洲精品粉嫩美女一区| 欧美+日韩+精品| 欧美在线一区亚洲| 精品国产超薄肉色丝袜足j| 成人av在线播放网站| av天堂在线播放| 久久精品国产99精品国产亚洲性色| 男女午夜视频在线观看| 法律面前人人平等表现在哪些方面| 久久精品国产清高在天天线| 一级毛片女人18水好多| 偷拍熟女少妇极品色| 亚洲国产精品合色在线| 男插女下体视频免费在线播放| 国内精品美女久久久久久| 亚洲人与动物交配视频| 精华霜和精华液先用哪个| 亚洲精品乱码久久久v下载方式 | 国产亚洲精品一区二区www| 欧美成人性av电影在线观看| 国产视频一区二区在线看| av女优亚洲男人天堂| 成人国产一区最新在线观看| 精品日产1卡2卡| 日韩欧美一区二区三区在线观看| 日韩大尺度精品在线看网址| 欧美日本视频| 国产一区二区在线av高清观看| 久久久久久久久中文| 51午夜福利影视在线观看| 精品福利观看| av女优亚洲男人天堂| 男女午夜视频在线观看| 欧美国产日韩亚洲一区| 99在线人妻在线中文字幕| 中文资源天堂在线| 国产精品久久视频播放| 免费在线观看影片大全网站| 首页视频小说图片口味搜索| 免费av毛片视频| 日日夜夜操网爽| 日韩欧美免费精品| 国产视频一区二区在线看| 在线观看免费午夜福利视频| av国产免费在线观看| 国产综合懂色| 午夜福利免费观看在线| 又黄又爽又免费观看的视频| 久久午夜亚洲精品久久| 色在线成人网| 最近最新中文字幕大全电影3| 国产精品日韩av在线免费观看| 国产伦在线观看视频一区| 国产精品 国内视频| 国产免费一级a男人的天堂| 桃色一区二区三区在线观看| 国产欧美日韩一区二区三| 国产精品野战在线观看| 久久6这里有精品| 精品一区二区三区视频在线观看免费| 免费在线观看成人毛片| 国产蜜桃级精品一区二区三区| 脱女人内裤的视频| 日本a在线网址| 99视频精品全部免费 在线| 亚洲人与动物交配视频| avwww免费| 啦啦啦观看免费观看视频高清| 天天添夜夜摸| av天堂中文字幕网| 可以在线观看的亚洲视频| 丁香欧美五月| 夜夜躁狠狠躁天天躁| 偷拍熟女少妇极品色| 麻豆成人av在线观看| 亚洲成人久久爱视频| 国产不卡一卡二| 人妻丰满熟妇av一区二区三区| 国产精品精品国产色婷婷| 国产黄片美女视频| 成人av在线播放网站| 啦啦啦免费观看视频1| xxx96com| 观看美女的网站| 日韩中文字幕欧美一区二区| 夜夜夜夜夜久久久久| 黄色女人牲交| 亚洲国产欧美网| 国产精品久久久久久人妻精品电影| 国产免费一级a男人的天堂| 少妇的丰满在线观看| 在线免费观看的www视频| 18+在线观看网站| 亚洲精华国产精华精| 久久精品人妻少妇| 蜜桃久久精品国产亚洲av| 亚洲精品国产精品久久久不卡| 又黄又爽又免费观看的视频| 在线观看日韩欧美| av欧美777| 免费一级毛片在线播放高清视频| 夜夜爽天天搞| 男人舔奶头视频| 丰满人妻熟妇乱又伦精品不卡| 丰满乱子伦码专区| 在线a可以看的网站| 美女 人体艺术 gogo| 男女做爰动态图高潮gif福利片| 国产av麻豆久久久久久久| 亚洲国产精品成人综合色| 欧美一区二区精品小视频在线| 精品99又大又爽又粗少妇毛片 | 国产私拍福利视频在线观看| 国产熟女xx| 免费观看的影片在线观看| 国产真实伦视频高清在线观看 | 亚洲第一电影网av| 久久久久久大精品| 99视频精品全部免费 在线| 桃色一区二区三区在线观看| 亚洲av美国av| 久久精品国产亚洲av涩爱 | 免费人成在线观看视频色| а√天堂www在线а√下载| 人妻夜夜爽99麻豆av| 成人亚洲精品av一区二区| 国产一区二区亚洲精品在线观看| 国产精品1区2区在线观看.| 亚洲人成伊人成综合网2020| 国产精品,欧美在线| 熟妇人妻久久中文字幕3abv| 午夜精品在线福利| 久久精品人妻少妇| 国产成人欧美在线观看| 国产一区二区在线观看日韩 | www日本在线高清视频| 超碰av人人做人人爽久久 | 国产成人系列免费观看| 精品99又大又爽又粗少妇毛片 | 两性午夜刺激爽爽歪歪视频在线观看| 欧美成人性av电影在线观看| 亚洲精品一卡2卡三卡4卡5卡| 中国美女看黄片| 老熟妇仑乱视频hdxx| 久久久国产成人精品二区| 欧美bdsm另类| 最近最新中文字幕大全电影3| 老汉色av国产亚洲站长工具| 人人妻人人看人人澡| 免费观看人在逋| 精品福利观看| 日韩亚洲欧美综合| 日本五十路高清| 国产成人啪精品午夜网站| 男女床上黄色一级片免费看| 中文字幕高清在线视频| 国产黄色小视频在线观看| 亚洲专区国产一区二区| 麻豆国产av国片精品| 亚洲av成人不卡在线观看播放网| 国产69精品久久久久777片| 亚洲片人在线观看| 国产成人a区在线观看| 国产一区二区三区在线臀色熟女| 99久久九九国产精品国产免费| 亚洲成av人片免费观看| 此物有八面人人有两片| 午夜福利成人在线免费观看| 日本黄色视频三级网站网址| 久久久国产精品麻豆| 欧美日韩国产亚洲二区| 午夜老司机福利剧场| 国产乱人视频| 国产成人欧美在线观看| 亚洲激情在线av| 国产色婷婷99| 欧美性猛交黑人性爽| 又黄又粗又硬又大视频| 嫩草影院入口| 每晚都被弄得嗷嗷叫到高潮| 国产日本99.免费观看| 亚洲美女视频黄频| 免费在线观看日本一区| 免费大片18禁| 一级黄色大片毛片| 精品熟女少妇八av免费久了| 一级a爱片免费观看的视频| 黄色丝袜av网址大全| 国产一区二区亚洲精品在线观看|