• <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
    免费在线观看视频国产中文字幕亚洲| 免费看日本二区| 中文字幕另类日韩欧美亚洲嫩草| 老汉色av国产亚洲站长工具| 亚洲av成人av| 视频区欧美日本亚洲| 黄色视频不卡| 免费高清视频大片| 哪里可以看免费的av片| 国产欧美日韩一区二区三| 宅男免费午夜| 男女床上黄色一级片免费看| 亚洲男人天堂网一区| 长腿黑丝高跟| 国产精品爽爽va在线观看网站 | 亚洲人成伊人成综合网2020| 日韩一卡2卡3卡4卡2021年| 黄色片一级片一级黄色片| 国产日本99.免费观看| 欧美丝袜亚洲另类 | 亚洲国产精品合色在线| 国产视频一区二区在线看| 少妇粗大呻吟视频| 97碰自拍视频| 18禁黄网站禁片免费观看直播| 国产精品久久久人人做人人爽| 成人欧美大片| 久久久久免费精品人妻一区二区 | 国产免费av片在线观看野外av| 波多野结衣高清作品| 啦啦啦 在线观看视频| 成人亚洲精品一区在线观看| 欧美性猛交╳xxx乱大交人| 亚洲七黄色美女视频| 十八禁人妻一区二区| 在线观看免费午夜福利视频| 亚洲欧美激情综合另类| xxxwww97欧美| 精品久久久久久成人av| 制服人妻中文乱码| 窝窝影院91人妻| 亚洲专区字幕在线| 黄色成人免费大全| 老司机靠b影院| 亚洲欧美日韩无卡精品| 日韩 欧美 亚洲 中文字幕| 久久精品国产清高在天天线| a级毛片a级免费在线| 最近最新免费中文字幕在线| 国产亚洲av嫩草精品影院| 亚洲成人久久爱视频| 特大巨黑吊av在线直播 | 岛国视频午夜一区免费看| 露出奶头的视频| 国产精品久久视频播放| 久久久国产成人免费| 亚洲第一欧美日韩一区二区三区| 久久午夜综合久久蜜桃| 精品国内亚洲2022精品成人| 美女国产高潮福利片在线看| 亚洲精品国产区一区二| 精品人妻1区二区| 久久精品夜夜夜夜夜久久蜜豆 | 国产精品久久久久久亚洲av鲁大| 国产在线观看jvid| 欧美午夜高清在线| 国产一区二区激情短视频| 免费高清在线观看日韩| 自线自在国产av| 亚洲美女黄片视频| 免费电影在线观看免费观看| 极品教师在线免费播放| 精品免费久久久久久久清纯| 久久亚洲真实| 给我免费播放毛片高清在线观看| 在线观看日韩欧美| x7x7x7水蜜桃| 精品日产1卡2卡| 欧美一级a爱片免费观看看 | 免费在线观看视频国产中文字幕亚洲| 亚洲熟女毛片儿| 午夜福利在线在线| 亚洲精品久久成人aⅴ小说| 亚洲熟妇中文字幕五十中出| 欧美日韩乱码在线| 国产成人影院久久av| 久久久精品国产亚洲av高清涩受| 亚洲精品国产一区二区精华液| 中文字幕久久专区| 亚洲国产日韩欧美精品在线观看 | 亚洲美女黄片视频| 久久亚洲精品不卡| 国产欧美日韩精品亚洲av| 国产高清有码在线观看视频 | 欧美zozozo另类| 午夜精品在线福利| 欧美丝袜亚洲另类 | 老鸭窝网址在线观看| 一个人免费在线观看的高清视频| 国产精品久久久久久人妻精品电影| 俺也久久电影网| 免费在线观看成人毛片| www国产在线视频色| 国产人伦9x9x在线观看| 人人澡人人妻人| 精品国内亚洲2022精品成人| 亚洲av成人一区二区三| 777久久人妻少妇嫩草av网站| 精品国产亚洲在线| 国产av在哪里看| 他把我摸到了高潮在线观看| av有码第一页| 一个人观看的视频www高清免费观看 | 桃红色精品国产亚洲av| 黄色a级毛片大全视频| 国产熟女xx| 不卡一级毛片| 欧美黑人精品巨大| 欧美日韩中文字幕国产精品一区二区三区| 免费搜索国产男女视频| 无人区码免费观看不卡| 亚洲第一电影网av| netflix在线观看网站| 亚洲中文字幕一区二区三区有码在线看 | 美女高潮喷水抽搐中文字幕| 午夜亚洲福利在线播放| 欧美色欧美亚洲另类二区| 亚洲中文日韩欧美视频| 老汉色∧v一级毛片| 久久精品国产亚洲av香蕉五月| 亚洲性夜色夜夜综合| tocl精华| www.精华液| 韩国av一区二区三区四区| 久久婷婷成人综合色麻豆| 在线免费观看的www视频| 老鸭窝网址在线观看| 99久久综合精品五月天人人| 国产成人精品久久二区二区免费| 国语自产精品视频在线第100页| 国产野战对白在线观看| 长腿黑丝高跟| 久久精品91无色码中文字幕| 日本撒尿小便嘘嘘汇集6| 亚洲欧洲精品一区二区精品久久久| 国产熟女xx| 老司机靠b影院| 欧美黑人精品巨大| 天天添夜夜摸| 人人妻人人澡人人看| 久久久久久久久免费视频了| 曰老女人黄片| 国产精品久久久久久亚洲av鲁大| 天堂√8在线中文| 这个男人来自地球电影免费观看| 成熟少妇高潮喷水视频| 老司机深夜福利视频在线观看| 亚洲 国产 在线| 男女床上黄色一级片免费看| 久久亚洲真实| 可以免费在线观看a视频的电影网站| 丁香六月欧美| 在线观看66精品国产| 天天躁夜夜躁狠狠躁躁| 18禁国产床啪视频网站| 国产成人精品久久二区二区91| 国产日本99.免费观看| 国产又黄又爽又无遮挡在线| 美国免费a级毛片| 国产爱豆传媒在线观看 | 国产一卡二卡三卡精品| 国产精品久久久久久精品电影 | 欧美最黄视频在线播放免费| 婷婷亚洲欧美| 亚洲专区中文字幕在线| 黄片小视频在线播放| 免费人成视频x8x8入口观看| 欧美丝袜亚洲另类 | 波多野结衣高清无吗| 国产精品,欧美在线| 两个人视频免费观看高清| av片东京热男人的天堂| 怎么达到女性高潮| 久久国产精品影院| 男人舔女人下体高潮全视频| 国产精品av久久久久免费| 亚洲人成电影免费在线| 人人澡人人妻人| 一区二区三区激情视频| av在线播放免费不卡| 狠狠狠狠99中文字幕| 夜夜夜夜夜久久久久| 人人妻人人看人人澡| 亚洲精品在线观看二区| 91成年电影在线观看| 欧美色视频一区免费| 日韩欧美国产在线观看| 久久久久国产精品人妻aⅴ院| 亚洲美女黄片视频| 丝袜美腿诱惑在线| 免费在线观看亚洲国产| 亚洲中文日韩欧美视频| 999久久久国产精品视频| 国产成年人精品一区二区| 午夜免费激情av| 午夜久久久久精精品| 久久国产精品影院| 日韩欧美一区二区三区在线观看| 日本a在线网址| 女人被狂操c到高潮| 91九色精品人成在线观看| 日韩免费av在线播放| 久久中文看片网| 18美女黄网站色大片免费观看| 亚洲男人的天堂狠狠| 国产免费男女视频| 亚洲美女黄片视频| 中文字幕另类日韩欧美亚洲嫩草| 亚洲欧美激情综合另类| 国产亚洲精品一区二区www| 欧美又色又爽又黄视频| 国产三级黄色录像| 国产成人av教育| 国产蜜桃级精品一区二区三区| 免费搜索国产男女视频| 法律面前人人平等表现在哪些方面| 少妇 在线观看| 此物有八面人人有两片| 国产精品久久久久久亚洲av鲁大| 91成年电影在线观看| 一级作爱视频免费观看| 日日摸夜夜添夜夜添小说| 18禁观看日本| 国产午夜精品久久久久久| 午夜精品久久久久久毛片777| 高潮久久久久久久久久久不卡| 动漫黄色视频在线观看| 黄色视频不卡| 别揉我奶头~嗯~啊~动态视频| 成人精品一区二区免费| 亚洲 国产 在线| 在线天堂中文资源库| 亚洲人成电影免费在线| 一级毛片高清免费大全| 波多野结衣高清无吗| 天天添夜夜摸| 国产97色在线日韩免费| 免费在线观看视频国产中文字幕亚洲| 亚洲成人精品中文字幕电影| 中文字幕另类日韩欧美亚洲嫩草| 精品久久蜜臀av无| www日本在线高清视频| 久久精品夜夜夜夜夜久久蜜豆 | 亚洲国产精品999在线| 两个人视频免费观看高清| 国产精品久久电影中文字幕| 午夜两性在线视频| 在线天堂中文资源库| 极品教师在线免费播放| 亚洲av电影不卡..在线观看| 少妇 在线观看| 亚洲男人天堂网一区| 国产真实乱freesex| netflix在线观看网站| 女人高潮潮喷娇喘18禁视频| 精品国产乱码久久久久久男人| 国产久久久一区二区三区| 麻豆成人av在线观看| 精品乱码久久久久久99久播| 禁无遮挡网站| 免费一级毛片在线播放高清视频| 日韩欧美国产一区二区入口| 91av网站免费观看| 色综合欧美亚洲国产小说| 少妇熟女aⅴ在线视频| 国产又色又爽无遮挡免费看| 久久精品91无色码中文字幕| 夜夜躁狠狠躁天天躁| 久久久久亚洲av毛片大全| 露出奶头的视频| 女人被狂操c到高潮| 国产乱人伦免费视频| 中文字幕人妻丝袜一区二区| 国产精品美女特级片免费视频播放器 | 欧美一区二区精品小视频在线| 国产精品日韩av在线免费观看| 亚洲人成伊人成综合网2020| 一本综合久久免费| 日韩欧美一区视频在线观看| 午夜免费激情av| 亚洲精品av麻豆狂野| 国产精品影院久久| 日韩欧美免费精品| 亚洲欧美精品综合久久99| 国产精品一区二区精品视频观看| 宅男免费午夜| 99热6这里只有精品| 手机成人av网站| 国产精品二区激情视频| 国产亚洲欧美98| 可以在线观看毛片的网站| 精品久久久久久,| 欧美黑人巨大hd| 欧美一级a爱片免费观看看 | avwww免费| 91国产中文字幕| 无限看片的www在线观看| 国产91精品成人一区二区三区| 国产国语露脸激情在线看| 日韩精品中文字幕看吧| 国产精品久久久久久精品电影 | 天天躁狠狠躁夜夜躁狠狠躁| www.www免费av| 中文字幕高清在线视频| 国产主播在线观看一区二区| 欧美色视频一区免费| 美国免费a级毛片| 亚洲成人国产一区在线观看| bbb黄色大片| www日本在线高清视频| 99在线人妻在线中文字幕| 亚洲欧美精品综合久久99| 亚洲一区二区三区色噜噜| 黄色视频,在线免费观看| 亚洲国产中文字幕在线视频| 亚洲aⅴ乱码一区二区在线播放 | √禁漫天堂资源中文www| 91国产中文字幕| 国产蜜桃级精品一区二区三区| 99热这里只有精品一区 | 国产高清视频在线播放一区| 88av欧美| 午夜福利视频1000在线观看| 亚洲成国产人片在线观看| 99精品在免费线老司机午夜| 波多野结衣高清无吗| 亚洲av日韩精品久久久久久密| 欧美性猛交╳xxx乱大交人| 夜夜爽天天搞| 久久草成人影院| 国产精品1区2区在线观看.| 岛国视频午夜一区免费看| 老鸭窝网址在线观看| 久99久视频精品免费| 亚洲第一欧美日韩一区二区三区| 免费一级毛片在线播放高清视频| 久久久久久九九精品二区国产 | 欧美激情极品国产一区二区三区| 国产国语露脸激情在线看| 国产激情偷乱视频一区二区| 国产精品综合久久久久久久免费| a级毛片a级免费在线| 91麻豆av在线| 亚洲专区国产一区二区| 男女床上黄色一级片免费看| 非洲黑人性xxxx精品又粗又长| 亚洲国产欧美网| 国产一区在线观看成人免费| 1024视频免费在线观看| 怎么达到女性高潮| 人人妻,人人澡人人爽秒播| 亚洲av熟女| 国产精品一区二区免费欧美| 欧美+亚洲+日韩+国产| 一级片免费观看大全| 怎么达到女性高潮| 国产又爽黄色视频| 日韩有码中文字幕| 欧美国产精品va在线观看不卡| 精品国产一区二区三区四区第35| 国产精品精品国产色婷婷| 97人妻精品一区二区三区麻豆 | 日韩av在线大香蕉| 免费在线观看亚洲国产| 男男h啪啪无遮挡| 久久久久久九九精品二区国产 | 午夜影院日韩av| 中文字幕另类日韩欧美亚洲嫩草| 中文字幕精品亚洲无线码一区 | 不卡一级毛片| av有码第一页| 国产成人影院久久av| 99国产综合亚洲精品| 一区二区三区精品91| 国产亚洲欧美精品永久| 亚洲va日本ⅴa欧美va伊人久久| 在线永久观看黄色视频| 日韩有码中文字幕| 曰老女人黄片| 变态另类丝袜制服| 欧美丝袜亚洲另类 | 中文字幕人成人乱码亚洲影| 久久九九热精品免费| 国产熟女xx| 黄片小视频在线播放| 久久久久久亚洲精品国产蜜桃av| 日本成人三级电影网站| 91成年电影在线观看| 校园春色视频在线观看| 身体一侧抽搐| 国产成人精品无人区| 日日摸夜夜添夜夜添小说| 日韩精品免费视频一区二区三区| 好男人在线观看高清免费视频 | 国产精品久久久久久人妻精品电影| 人妻丰满熟妇av一区二区三区| 少妇粗大呻吟视频| 久久这里只有精品19| 国产精品久久久久久亚洲av鲁大| 欧美丝袜亚洲另类 | 精品人妻1区二区| 99久久精品国产亚洲精品| 大型黄色视频在线免费观看| 一个人观看的视频www高清免费观看 | 夜夜躁狠狠躁天天躁| 久热爱精品视频在线9| 他把我摸到了高潮在线观看| av超薄肉色丝袜交足视频| 成人18禁在线播放| 精品少妇一区二区三区视频日本电影| 在线视频色国产色| 亚洲国产欧美一区二区综合| 国产精品久久电影中文字幕| 国产视频一区二区在线看| 丝袜美腿诱惑在线| 人人澡人人妻人| 国产又爽黄色视频| 国产成人av教育| 亚洲熟女毛片儿| 91av网站免费观看| 看免费av毛片| 亚洲精品久久成人aⅴ小说| 非洲黑人性xxxx精品又粗又长| 精品不卡国产一区二区三区| 久久人妻福利社区极品人妻图片| xxx96com| 在线观看一区二区三区| 美女 人体艺术 gogo| 国产精品乱码一区二三区的特点| 色综合站精品国产| 制服人妻中文乱码| 亚洲真实伦在线观看| 亚洲五月婷婷丁香| 欧美性长视频在线观看| 变态另类丝袜制服| 精品久久久久久久久久久久久 | 免费女性裸体啪啪无遮挡网站| 国产1区2区3区精品| 18禁国产床啪视频网站| 午夜成年电影在线免费观看| 亚洲av成人不卡在线观看播放网| 操出白浆在线播放| 亚洲美女黄片视频| 国产片内射在线| 老汉色av国产亚洲站长工具| 国产亚洲精品综合一区在线观看 | 色av中文字幕| 1024香蕉在线观看| 麻豆一二三区av精品| 免费女性裸体啪啪无遮挡网站| 久久精品国产综合久久久| 亚洲久久久国产精品| 日本三级黄在线观看| www.自偷自拍.com| 搡老妇女老女人老熟妇| 精品一区二区三区av网在线观看| 国产野战对白在线观看| 夜夜爽天天搞| 日日夜夜操网爽| 村上凉子中文字幕在线| 国产成人精品无人区| 怎么达到女性高潮| 村上凉子中文字幕在线| 两个人视频免费观看高清| 国产亚洲精品av在线| 亚洲国产看品久久| 老司机午夜十八禁免费视频| 国产高清视频在线播放一区| 亚洲五月天丁香| 又黄又粗又硬又大视频| 成人永久免费在线观看视频| 成人三级黄色视频| 中国美女看黄片| 观看免费一级毛片| 婷婷精品国产亚洲av在线| 一边摸一边抽搐一进一小说| 亚洲精品av麻豆狂野| 美女大奶头视频| 午夜福利视频1000在线观看| av电影中文网址| 久久久久久久久久黄片| 97超级碰碰碰精品色视频在线观看| 亚洲欧美精品综合一区二区三区| 国产精品爽爽va在线观看网站 | 啦啦啦观看免费观看视频高清| 免费看a级黄色片| 满18在线观看网站| 制服人妻中文乱码| 国产黄片美女视频| 国产一卡二卡三卡精品| aaaaa片日本免费| 又紧又爽又黄一区二区| 午夜亚洲福利在线播放| 午夜福利成人在线免费观看| 精品国产国语对白av| 97人妻精品一区二区三区麻豆 | 中国美女看黄片| 亚洲五月婷婷丁香| 久99久视频精品免费| e午夜精品久久久久久久| svipshipincom国产片| 天堂影院成人在线观看| 国产精品98久久久久久宅男小说| 免费看a级黄色片| 欧美激情 高清一区二区三区| 狠狠狠狠99中文字幕| 国产极品粉嫩免费观看在线| 亚洲国产中文字幕在线视频| 男女那种视频在线观看| netflix在线观看网站| 日韩精品中文字幕看吧| 国产男靠女视频免费网站| 少妇被粗大的猛进出69影院| 岛国视频午夜一区免费看| 日本免费a在线| 禁无遮挡网站| 99久久国产精品久久久| 欧美黄色淫秽网站| 在线播放国产精品三级| 欧美日本视频| 大型黄色视频在线免费观看| 一区二区日韩欧美中文字幕| 美女高潮到喷水免费观看| 在线观看免费午夜福利视频| 国产单亲对白刺激| 黄片播放在线免费| 久久精品亚洲精品国产色婷小说| 精品国产超薄肉色丝袜足j| 伊人久久大香线蕉亚洲五| 在线永久观看黄色视频| 国产单亲对白刺激| 亚洲真实伦在线观看| 在线观看日韩欧美| 日本黄色视频三级网站网址| 久久久久久久精品吃奶| 熟妇人妻久久中文字幕3abv| 亚洲av日韩精品久久久久久密| 夜夜躁狠狠躁天天躁| 亚洲第一青青草原| √禁漫天堂资源中文www| 搡老岳熟女国产| 欧美黄色淫秽网站| 国产av在哪里看| 18禁黄网站禁片午夜丰满| 在线观看www视频免费| cao死你这个sao货| 搞女人的毛片| 极品教师在线免费播放| 可以在线观看毛片的网站| 国产精品二区激情视频| 日韩欧美国产在线观看| 久久亚洲精品不卡| 中出人妻视频一区二区| 日韩欧美国产一区二区入口| 免费在线观看黄色视频的| 色哟哟哟哟哟哟| 午夜福利视频1000在线观看| 老司机靠b影院| 免费在线观看亚洲国产| 精华霜和精华液先用哪个| 法律面前人人平等表现在哪些方面| 亚洲片人在线观看| 国产精品永久免费网站| 国产精品,欧美在线| 黄片播放在线免费| 亚洲精品国产精品久久久不卡| 精品久久久久久久人妻蜜臀av| 国产一区二区激情短视频| 麻豆一二三区av精品| 色综合站精品国产| 午夜日韩欧美国产| 午夜福利免费观看在线| 日日爽夜夜爽网站| videosex国产| 在线观看舔阴道视频| 欧美激情高清一区二区三区| 中文资源天堂在线| 少妇裸体淫交视频免费看高清 | 亚洲成国产人片在线观看| 亚洲精品一区av在线观看| 成人国语在线视频| 午夜福利18| 亚洲成a人片在线一区二区| 波多野结衣高清无吗| x7x7x7水蜜桃| 国产av在哪里看| xxxwww97欧美| 免费看a级黄色片| 少妇的丰满在线观看| 搡老岳熟女国产| 亚洲avbb在线观看| 日韩欧美一区视频在线观看| 在线播放国产精品三级| 深夜精品福利| 久久国产亚洲av麻豆专区| 成在线人永久免费视频| 欧美乱色亚洲激情| 欧美日韩乱码在线| 国产精品98久久久久久宅男小说| 十八禁人妻一区二区| 悠悠久久av| 很黄的视频免费| 97人妻精品一区二区三区麻豆 | 国产精品亚洲美女久久久| 国产成人精品无人区|