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

    Lattice Boltzmann Flux Solver:An Efficient Approach for Numerical Simulation of Fluid Flows*

    2014-04-24 10:53:08ShuChangWangYangWu

    Shu Chang,Wang Y,Yang L M,Wu J

    1.Department of Mechanical Engineering,National University of Singapore,10Kent Ridge Crescent,Singapore,119260;2.College of Aerospace Engineering,Nanjing University of Aeronautics and Astronautics,Nanjing,210016,P.R.China(Received 20January 2014;revisied 18February 2014;accepted 22February 2014)

    1 Introduction

    Computational fluid dynamics(CFD)is to apply a numerical method to solve governing equations of fluid flows on the computer.Among various numerical methods available[1-8],the finite volume method(FVM)is the most popular approach in CFD.This is because numerical discretization by FVM is in line with application of physical conservation laws to a control cell.The discrete forms of governing equations by FVM usually involve the conservative variables at cell centers and numerical fluxes at cell interfaces.From numerical point of view,only the conservative variables at cell centers are defined as unknowns,which can be given from the solution of discrete governing equations.In the solution process,we need to use conservative variables at cell centers to evaluate numerical fluxes at cell interfaces.This process is often termed flux solver.Currently,there are three major flux solvers in CFD.One is based on the smooth function approximation.In this solver,a smooth function,which could be a polynomial[9]or a radial basis function[10],is applied to approximate the solution in the local region.The coefficients in the smooth function can be determined by collocation method.Once the smooth function is decided,its integral or derivative can be given in a straightforward way.It should be noted that this solver is a mathematical approach,which can be applied to general engineering problems.However,this solver cannot resolve discontinuity problems such as compressible flows with shock wave.To re-solve shock wave problems in CFD,the Riemann solver or approximate Riemann solver is often used.The pioneer work in this category was made by Godunov[11],who simplified the compressible flow into a series of Riemann problems and then solved one-dimensional(1D)Euler equations to get local solution.After the work of Godunov[11],various approximate Riemann solvers were presented[12-16].These solvers usually pursue approximate solution of 1DEuler equations along the normal direction to the cell interface.Thus,they can only be used to evaluate inviscid flux.For compressible viscous flows,the viscous flux is still evaluated by the smooth function approximation.In the literature,there is another type of flux solver called gas kinetic flux solver[17-19],which evaluates inviscid and viscous fluxes simultaneously from local solution of multi-dimensional Boltzmann equation.The solvers in this category can be well applied to simulate both incompressible and compressible flows.But they are usually more complicated and less efficient than the smooth function-based solvers and Riemann solvers.In this paper,we will present a new flux solver,which is based on local solution of lattice Boltzmann equation(LBE).

    In recent years,lattice Boltzmann method(LBM)[20-32]has received more and more attention due to its simplicity,easy implementation and parallel nature.In LBM,the density distribution functions are taken as unknowns and LBE is an algebraic formulation.Once the density distribution functions are known at a physical location,the macroscopic flow variables such as density and velocity can be easily computed from local conservation laws of mass and momentum.No differential equation and solution of algebraic equations are involved in the LBE solver.On the other hand,it is indicated that LBE solvers also suffer from some drawbacks.Due to uniformity of the lattice,the standard LBE solver is limited to the simple geometry and uniform mesh.For complex geometry and application on the non-uniform mesh,additional efforts such as interpolation have to be incorporated.The process may in-crease the complexity of the solver,and requires additional computational effort and virtual storage.The second drawback is the tie-up of time interval with mesh spacing.This drawback makes the adaptive and multi-block computation of LBE solvers extremely complicated.In addition,LBE solvers need more memory to store density distribution functions than the Navier-Stokes(N-S)solvers.Another drawback is that LBE solvers can only be applied to simulate viscous flows.Furthermore,the physical boundary conditions such as given pressure cannot be implemented directly in the LBE solver.As will be shown in this paper,all the above drawbacks of LBE solvers are completely removed by the lattice Boltzmann flux solver(LBFS).

    LBFS is based on Chapman-Enskog(C-E)expansion analysis,which is a bridge to link N-S equations and LBE.Usually,the C-E analysis is applied in the whole flow domain to verify that the macroscopic flow variables obtained by LBE solvers at any physical location and any time level can satisfy N-S equations.On the other hand,it was found that the C-E analysis can be applied at any location within a small streaming step.This idea has been well applied by Xu[17]in the development of gas kinetic scheme,where the flux at the cell interface is computed by local solution of BGK equation.In this work,the numerical fluxes at the cell interface are evaluated by local reconstruction of LBE solution.Two versions of LBFS are presented in this work.One is to locally apply 1Dcompressible LB model along the normal direction of cell interface for evaluation of inviscid flux.This version is only applicable for simulation of compressible inviscid flows.The other is to locally apply multi-dimensional LB model at the cell interface for evaluation of viscous and inviscid fluxes simultaneously.The performance of present LBFS will be investigated through some test examples.Numerical results demonstrate that LBFS can accurately and effectively simulate fluid flows with curved boundary and non-uniform mesh.It also removes the drawbacks of conventional LBM.

    2 Lattice Boltzmann Flux Solver(LBFS)for Compressible Inviscid Flows

    The integral form of Euler equations without source term can be written as

    where the conservative flow variables Wand inviscid flux Fnare given by

    whereρand pare the density and pressure of the mean flow,respectively.U=(u,v,w)is the velocity vector in the Cartesian coordinate system and n=(nx,ny,nz)denotes the unit normal vector on the control surface.Unrepresents the normal velocity,which is defined as the scalar product of the velocity vector and the unit normal vector,i.e.

    Eis the total energy of the mean flow,which is defined as

    Here e=p/[(γ-1)ρ]is the potential energy of the mean flow,andγis the specific heat ratio.On the control surface,the tangential velocity Uτ= (Uτx,Uτy,Uτz)can be computed by

    Applying Eq.(1)to a control volume gives

    where Iis the index of a control volume,ΩIand Nfrepresent the volume and the number of the faces of the control volume I.dSidenotes the area of the ith face of the control volume.As indicated in the introduction,the flux solver needs to reconstruct numerical flux Fnat each cell interface from the conservative variables WIat cell centers.In this section,Fnwill be computed from the solution of 1Dcompressible LB model to a local Riemann problem.When 1DLB model is applied along the normal direction to the cell interface,only density,pressure and normal velocity are involved.Thus,before we address how to apply the 1Dcompressible LB model to reconstruct Fn,it is better to rewrite expression of Fnin terms of density,pressure,normal velocity and tangential velocity.From Eq.(5),we can express the velocity components in the Cartesian coordinate system in terms of normal velocity and tangential velocity as

    Using Eq.(7)and the expression of potential energy,Fncan be rewritten as

    It can be seen clearly from Eq.(8)that,to evaluate numerical flux Fn,we need to know the density,pressure,normal velocity and tangential velocity at the cell interface.This task can be fulfilled by local application of 1Dcompressible LB model to the Riemann problem defined at cell interface.In this work,the non-free parameter D1Q4LB model presented in Refs.[33-34]is adopted.This model is derived from conservation forms of moments,which can be used to simulate hypersonic flows with strong shock waves.The non-free parameter D1Q4model is shown in Fig.1.The equilibrium distribution functions and lattice velocities of this model are given below,where giis the equilibrium distribution function in the ith direction of phase space,diis the lattice velocity in the ith direction,c is the peculiar ve-locity of particles defined as c=(Dis the dimension of space).Note that when the above 1Dmodel is applied along the normal direction to the cell interface,u has to be replaced by Un.

    Fig.1 Configuration of non-free parameter D1Q4model

    Next,we will show how to apply the nonfree parameter D1Q4model to evaluate Fnat cell interface.As shown in Fig.1,at any physical location,D1Q4model has 4moving particles.Now,we consider a local Riemann problem around a cell interface as shown in Fig.2.To compute Fn,we need to know distribution functions of 4moving particles at the cell interface.In the framework of LBM,the moving particles are actually streamed from neighbouring points.As illustrated in Fig.3,by giving a streaming step δt,particles 1and 3from left side of interface will stream to the cell interface while particles 2and 4 from right side of interface will also stream to the cell interface.Mathematically,the streaming process provides the distribution functions of four moving particles at cell interface as

    Fig.2 Configuration of a Riemann problem

    Fig.3 Streaming process of D1Q4model at the cell interface

    where gLiand gRiare the equilibrium distribution functions at the left and right sides of cell interface.For the Riemann problem,they are given from information at left and right cell centers.With flow variables,they can be computed by using Eq.(9).With Eq.(11),there are two basic ways to evaluate the numerical flux Fnat the cell interface.The first way is to compute the flow variables(density,pressure and normal velocity)first,and then substitute them into Eq.(8)to compute Fn.The density,normal velocity and pressure can be computed by

    where eiis the lattice velocity,e1=d1,e2=-d1,e is the po-tential energy of particles(Dis the dimension of space and takes 1for the 1Dmodel).The tangential density Uτat the cell interface can be given from mean value of ULτand URτ,where ULτand URτare the tangential velocity at the left and right side of cell interface,respectively.Alternatively,it can be approximated by

    Once the density,pressure,normal velocity and tangential velocity at the cell interface are computed by Eqs.(12-13),they can be substituted into Eq.(8)to compute Fn.This way is equivalent to use equilibrium distribution functions at the cell interface to compute Fn.From CE analysis,this way has very little numerical dissipation,which may not be able to get stable solution for problems with strong shock waves.To compute Fnwith numerical dissipation,we can use distribution function given in Eq.(11)to compute Fndirectly.In fact,ρUnin Fnhas been calculated by Eq.(12b).Other terms in Fncan be computed by the following formulations

    Similar to Eq.(13),ρUnUτandρUn|Uτ|2can be approximated by

    Overall,the basic solution procedure of this LBFS can be summarized below:

    (1)At first,we need to choose a 1DLB model such as non-free parameter D1Q4model.The LB model provides expressions for equilibrium distribution functions and lattice velocities.

    (2)For the considered cell interface with unit normal vector n= (nx,ny,nz),obtain flow variables(density,pressure,velocity components)at the left and right sides of interface from two neighbouring cell centers(MUSCL interpolation with limiter may be used for high-order schemes).Then use Eqs.(3,5)to calculate the normal and tangential velocities at the left and right sides of interface.

    (3)Use Eq.(9)to calculate gL1,gL3,gR2,gR4by using density,pressure and normal velocity.

    (4)Compute the density,normal velocity,pressure and tangential velocity at the cell interface by using Eqs.(12-13),and then substitute them into Eq.(8)to calculate numerical flux Fn.Alternatively,use Eqs.(12b,14-17)to compute Fndirectly(this way is recommended for hypersonic flows with strong shock waves).

    (5)Once numerical fluxes at all cell interfaces are obtained,solve ordinary differential equations(6)by using 4-stage Runge-Kutta scheme.

    For simulation of viscous flows,one also needs to use a smooth function to approximate the viscous flux.

    3 Lattice Boltzmann Flux Solver(LBFS)for Incompressible Flows

    From C-E expansion analysis[22,32],the incompressible Navier-Stokes(N-S)equations

    can be recovered by the following LBE

    whereρis the fluid density,u the flow velocity and pthe pressure.r represents a physical location,τis the single relaxation parameter;fαis the density distribution function along theαdirection;feqαis its corresponding equilibrium state;δtis the streaming time step and eαis the particle velocity in theαdirection;Nis the number of discrete particle velocities.The relationships between the density distribution functions and flow variables as well as fluxes in the N-S equations are

    whereβandγrepresent the space coordinate directions,and eαβis the component of the lattice velocity vector eαin theβ-coordinate direction.As shown in Refs.[22,32],to recover N-S equations by Eq.(20),εf(1)αcan be approximated by

    Substituting Eq.(24)into Eq.(23)gives

    The equilibrium distribution functiondepends on the lattice velocity model used.For example,when the following two-dimensional D2Q9lattice velocity model

    where c=δx/δt,δxis the lattice spacing.For the case ofδx=δt,which is often used in the literature and also adopted in this work,c is taken as 1.The coefficients wαand the sound speed csare given as:w0=4/,w1=w2=w3=w4=1/9and w5=w6=w7=w8=1/36.cs=.The relaxa-tion parameterτis linked to the kinematic viscosity of fluid through C-E expansion analysis by the following relationship

    The pressure can be calculated from the equation of state by

    Using Eqs.(22)and(23),for the two-dimensional case,Eqs.(18)and(19)can be rewritten as

    where

    When a cell-centered FVM is applied to solve Eq.(30),the flow propertiesρandρu at the cell center can be obtained by marching in time.The fluxes at the cell interface can be evaluated by local reconstruction of LBM solution.By integrating Eq.(30)over a control volumeΩi,we have

    whereΔViis the volume ofΩi,andΔSkis the area of the kth control surface enclosingΩi.nxand nyare the xand ycomponents of the unit outward normal vector on the kth control surface.Obviously,once the fluxes at all cell interfaces are known,Eq.(34)can be solved by well established numerical schemes such as the 4-stage Runge-Kutta method.Thus,the evaluation of flux Rkat the cell interface is the key in the solution process.The detailed expression of Rkdepends on the lattice velocity model.By defining fas

    When the D2Q9lattice velocity model is used,Rkcan be written in detail as follows

    Obviously,the key issue in the evaluation of the flux Rkis to perform an accurate evaluation of fandat the cell interface.In the following,we will show the detailed calculation ofand fat the cell interface.

    Consider a cell interface between two control cellsΩiandΩi+1as shown in Fig.4.It is assumed that the physical location for the two cell centers and their interface is ri,ri+1and r respectively.

    Fig.4 Local reconstruction of LBM solution at a cell interface

    Using Taylor series expansion,we have

    From Eqs.(37)and(24),we can get the following form

    Eq.(38)shows that once we have the equilibrium distribution functions(r,t),(r-eαδt,t-δt)at the cell interface and its surrounding points,we can have the full information of distribution function at the interface.Note that the approximation for Eq.(38)is the second order of accuracy inδt.Using Eq.(27),the equilibrium distribution function feqαcan be computed from the fluid densityρand flow velocity u.With the given density and velocity at the cell center,the respective density and velocity at location(r-eαδt)can be easily obtained by interpolation.One of interpolation forms can be written as

    With computedρ(r-eαδt)and u(r-eαδt)by Eqs.(39,40),(r-eαδt,t-δt)can be calculated by Eq.(27).Now,we are only left to determine(r,t)as shown in Eq.(38).Again,with Eq.(27),the calculation of(r,t)is equivalent to computingρ(r,t)and u(r,t).Using Eqs.(21,22),the conservative variablesρandρu can be computed by

    Since fαcan be written as,application of Eq.(20)at the cell interface leads to

    Furthermore,by substituting Eq.(38)into Eq.(43),we obtain

    Equation(44)is actually equivalent to fα(r,t)=r,t)+(r,t).Finally,Summation of Eq.(44)overαand applying the compatibility condition gives

    Eqs.(45,46)show that the conservative flow variables at the cell interface are fully determined from the equilibrium distribution functions at the surrounding points.As equilibrium distribution functions only depend on the macroscopic flow variables,there is no need to store the densi-ty distribution functions for all the time levels.In fact,at any time step,we locally reconstruct a LBM solution at each cell interface independently.The reconstruction process is applied locally and repeated from one time level to another time level.Overall,the basic solution procedure of LBFS can be summarized below:

    (1)At beginning,we need to choose a lattice velocity model such as D2Q9model.Then we need to specify a streaming time stepδt.The choice ofδtshould satisfy the constraint that the location of(r-eαδt)must be within either the cell Ωior the cellΩi+1.Note that as local LBM solution is reconstructed at each cell interface,different interfaces could use differentδt.This provides agreat flexibility for application if we use non-uniform mesh or solve problems with a curved boundary.Onceδtis chosen,the single relaxation parameterτin LBFS is calculated by Eq.(28).

    (2)For the considered interface position r,identify its surrounding positions(r-eαδt),and then use Eqs.(39,40)to compute the macroscopic flow variables at those positions.

    (3)Use Eq.(27)to calculate the equilibrium density distribution function(r-eαδt,t-δt).

    (4)Compute the macroscopic flow variables at the cell interface by using Eqs.(45)and(46),and further calculate(r,t)by Eq.(27).

    (6)Compute the fluxes at the cell interface by Eq.(36).

    (7)Once fluxes at all cell interfaces are obtained,solve ordinary differential Eq.(34)by using 4-stage Runge-Kutta scheme.

    It is indicated that the present LBFS can be used to simulate both incompressible viscous flows and incompressible inviscid flows.For the inviscid flow,we just simply setτ=0.5.Another point to note is that the time marching step used in solving Eq.(34)and the streaming time stepδtused in LBFS are independent.δtcan be selected differently at different interface and dif-ferent time level.Numerical experiments show thatδthas no effect on the solution accuracy.

    4 Numerical Examples and Discussion

    In this section,the developed LBFS is validated by its application to solve some test problems.In all following simulations,the non-free parameter D1Q4model[33-34]is used for simulation of compressible inviscid flows,and the D2Q9lattice velocity model is applied for simulation of two-dimensional incompressible viscous flows.

    4.1 Simulation of two-dimensional compressible inviscid flows

    At first,the LBFS developed in Section 2 will be applied to simulate three two-dimensional compressible inviscid flows.They are the flow around a NACA0012airfoil,the flow around a forward facing step,and the flow around a circular cylinder.For the flow around the NACA0012airfoil,the free-stream Mach number is taken as 0.8 and the angle of attack is chosen as 1.25°.Unstructured grid with 10 382cells is used for numerical computation.Both LBFS and Roe scheme are applied to solve this problem on the same computational mesh.It was found that the pressure coefficient distributions obtained by LBFS and Roe scheme are close to each other.The lift and drag coefficients(Cland Cd)obtained by LBFS are respectively 0.304 1and 0.023 7,which agree well with the results given from Roe scheme(Cl=0.283 6,Cd=0.021 5)and those of Stolcis and Johnston[35](Cl=0.339 7,Cd=0.022 8).Fig.5shows the pressure contours around the airfoil.As can be seen clearly,the shock wave on the upper surface is well captured by present solver.The second test example in this part is a stationary flow(Mach number equals 3)hitting a rectangular step.This problem has been well studied by Woodward and Colella[36],and is often used to investigate performance of new numerical methods for capturing the shock waves.In our computation,a uniform mesh size of 300×100is used.Fig.6shows the density contours computed by present solver.Our results are in good agreement with those in Ref.[36].It is noted that no special treatment around step corner is made in the present computation,which is often needed by conventional schemes.To further explore the capability of present solver for simulation of hypersonic flows with strong shock waves,the flow around a circular cylinder is simulated.For this case,a uniform mesh size of 160×40in the cylindrical coordinate system is used.It is well known that for this problem,conventional numerical schemes such as Roe scheme may encounter the″carbuncle phenomenon″in front of cylinder when the free stream Mach number is high.The″carbuncle phenomenon″may be due to unsatisfying of entropy condition and negative value of density in the local region.We have used different free-stream Mach numbers to test simulation of this problem by LBFS.For all the cases tested(free-stream Mach number up to 100),no″carbuncle phenom-enon″was found in the present results.This can be seen clearly from Fig.7,which shows pressure contours of Ma=3and 100.Both results show regular pressure distribution around the cylinder.

    Fig.5 Pressure contours around NACA0012airfoil

    Fig.6 Density contours for flow around a forward facing step

    Fig.7 Pressure contours for flow around a circular cylinder

    4.2 Simulation of compressible inviscid flows around ONERA M6wing

    Fig.8 Partial view of computational mesh for flow around ONERA M6wing

    To investigate the capability of present LBFS for solving practical flow problems,the threedimensional(3D)transonic flow around the ONERA M6wing is simulated.This is also a standard test case for 3Dcomputations.For numerical simulation,the free-stream Mach number is taken as 0.839 5and the angle of attack is chosen as 3.06°.The part of computational mesh is shown in Fig.8,which has 294 912cells.The pressure contours obtained by present solver are displayed in Fig.9.The″λ″shape shock wave on the upper surface of the wing can be seen clearly in Fig.9,which is in line with the result in Ref.[37].The pressure coefficient distribution at a section of z/b=0.65is shown in Fig.10.Also included in Fig.10are the experimental data given in Ref.[38].As can been clearly,the present results quantitatively compare very well with the experimental data.

    Fig.9 Pressure contours around ONERA M6wing

    Fig.10 Pressure coefficient distribution at section of z/b=0.65on M6wing

    4.3 Simulation of incompressible lid-driven flow in a square cavity

    The lid-driven flow in a square cavity is a standard test case for validating new numerical methods in simulation of incompressible viscous flows.The flow pattern of this problem is governed by the Reynolds number defined by Re=UL/ν,where Uis the lid speed,Lis the length of the cavity,andνis the kinematic viscosity of fluid.Two cases of this problem at moderate and high Reynolds numbers of 3 200and 7 500are considered in this work.LBFS introduced in Section 3will be applied to solve this problem and the following problems.

    To conduct numerical simulations,the nonuniform grid is generated according to the following formulation

    where Nand Mare the total number of mesh points in the xand ydirections respectively.With Eq.(47),the non-uniform grids of 101×101for Re=3 200and 121×121for Re=7 500are used respectively.In the present study,we set U=0.1 and L=1.The initial flow field is at rest.

    Table 1compares the locations of the primary vortex centers at Re=3 200and 7 500obtained by LBFS with those given by Ghia et al[39].As can be seen,the maximum relative error between present results and those of Ghia et al[39]is less than 1.1%.Fig.11displays u-velocity profile along the horizontal centerline and v-velocity profile along the vertical centerline of the considered two cases.As can be seen from this figure,the present results agree very well with those of Ghia et al[39].Fig.12shows the streamlines of Re=3 200,7 500.The most striking aspect of this figure is that the Reynolds number apparently has unique effect on flow patterns.Secondary and tertiary vortices appear and evolve into larger ones as Re becomes large.These results and observations are in good agreement with those of Ghia et al[39].

    Table 1 Locations of primary vortex centers at different Reynolds numbers

    Fig.11 uand vvelocity profiles along horizontal and vertical centerlines for a lid-driven cavity flow at Re=3 200,7 500

    Fig.12 Streamlines of a lid-driven cavity flow at Re=3 200,7 500

    Note that for this test example,we have also studied the effect of streaming distance in local reconstruction of LBM solution.It was found that when the streaming distance is less than half of mesh spacing in the two neighboring cells(this constraint guarantees that only interpolation is performed in each cell),any value of streaming distance will have no effect on the accuracy of solution.This is an appealing feature,which ensures that LBFS can be easily applied on non-uniform mesh.

    4.4 Simulation of incompressible polar cavity flow

    Although the complex lid-driven cavity flows have been successfully simulated to validate the present solver,the geometry of the cavity which only involves straight boundaries is nevertheless simple.To further illustrate the capability of LBFS for problems with curved boundary,apolar cavity flow is simulated on body-fitted meshes.The schematic diagram and the typical non-uniform mesh for this problem are depicted in Fig.13.As shown in Fig.13,a sector with an angle ofθ=1is bounded by two straight walls and two curved walls with radii of Riand Ro.The inner curved wall rotates with an azimuthal velocity of Uθ.The flow pattern of this problem is governed by the Reynolds number defined as Re=UθRi/ν.In this study,two cases of Re=60and 1 000are considered,and the following parameters are applied:Ri=1.0,Ro=2.0,ρ0=1.0and Uθ=0.1.Initially,the flow field is at rest.Fig.14shows the radial(ur)and azimuthal(uθ)velocity profiles along the horizontal line ofθ=0.5at Re=60and 1 000.The experimental results of Fuchs and Tillmark[40]and the numerical solutions of Shu et al[41]obtained by applying Taylor series expansion-and least-square-based LBM(TLLBM)are also included for comparison.Note that the present results and TLLBM results are both obtained on the same non-uniform grids,i.e.,61×61for Re=60and 81×81for Re=1 000.It can be seen that good agreements have been achieved between the present results and those of Fuchs and Tillmark[40]and Shu et al[41],which validate the reliability of the present solver for problems with curved boundary and use of non-uniform grid.The streamlines are shown in Fig.15.As can be seen,with increase of the Reynolds number,the primary vortex moves upward and reduces its size.At the same time,the two secondary vortices at the upper-right and lower-right corners enlarge their size.These observations agree well with those of Fuchs and Tillmark[40].

    Fig.13 Schematic diagram and a typical body-fitted mesh for flow in a polar lid-driven cavity

    4.5 Simulation of flow induced by an impulsively started cylinder

    In this part,LBFS is applied to simulate the unsteady flow induced by an impulsively started circular cylinder.The Reynolds number of this flow is defined as Re=UD/ν,where Uis the freestream velocity and Dis the diameter of cylinder.A wide range of Reynolds numbers from 102to 104are considered in this study to further demonstrate the capability of LBFS for effective simula-tion of unsteady flows at high Reynolds numbers.

    Fig.14 Comparison of radial(ur)and azimuthal(uθ)velocity profiles along the horizontal line of θ=0.5for the polar cavity flow

    Fig.15 Streamlines for the polar cavity flow at Re=60and 1 000

    In the present simulation,for flows at Re=550and 3 000,a mesh size of 301×201is used and the outer boundary is placed at 15diam-eters away from the cylinder center.For the flow at Re=9 950,the computational mesh is set as 301×351and the outer boundary is set as 4diameters away from the cylinder center.The flow parameters are set as:ρ=1.0,U=0.1and a=0.5.Initially,the flow field is at rest.

    For incompressible flows around the circular cylinder at high Reynolds number,apair of primary symmetric vortices will be developed at the rear of cylinder initially.With increase of the Reynolds number,the size of the two vortices is decreased.As time increases,the primary vortices will move away and detach from the rear of cylinder.In the meantime,apair of secondary symmetric vortices appears and becomes larger and stronger.The vortex structures exhibit the so-called″α″and″β″patterns.All these features have been well captured in present simulation.To save the space,these results are not displayed in this paper.Fig.16shows a quantitative comparison of the time evolution of the vortex length with experimental data of Bouard and Coutanceau[42].Obviously,good agreement has been achieved.For Re=550,the vortex length almost grows linearly with respect to time.For high Reynolds numbers(Re=3 000and 9 500),a slow increase in vortex length,which corresponds to the″fore-wake″region,can be observed when t<3.0s.When t>3s,a fast growth of the vortex length can be seen due to destruction of the″fore-wake″.Fig.17further compares the radial velocity along the symmetric axis at Re=3 000with experimental data of Bouard and Coutanceau[42]and numerical results of Niu et al[43].Once again,good agreement is achieved.

    Fig.16 Comparison of the vortex length for flow induced by impulsively started cylinder at different Reynolds numbers

    Fig.17 Comparison of the radial velocity along symmetric axis for flow induced by impulsively started cylinder at Re=3 000

    5 Conclusions

    In this paper,the LBFS is presented for simulation of compressible and incompressible flows.The solver is based on numerical discretization of FVM to the governing differential equations(Navier-Stokes equations or Euler equations).Specifically,the conservative flow variables at cell centers are given from the solution of discrete governing equations but numerical fluxes at cell interfaces are evaluated by local reconstruction of LBE solution from flow variables at cell centers.Two versions of LBFS are presented in this paper.One is to locally apply 1Dcompressible LB model along the normal direction to the cell interface for simulation of compressible inviscid flows.The other is to locally apply incompressible LB model at the cell interface for simulation of incompressible viscous flows.

    The present LBFS is well validated by its application to simulate some two-and three-dimensional compressible inviscid flows,and twodimensional incompressible viscous flows.Numerical results show that the compressible version of LBFS can well simulate compressible inviscid flows with strong shock waves,and its incom-pressible version can accurately simulate incompressible viscous flows with curved boundary and non-uniform mesh.It removes the drawbacks of conventional LBM such as limitation to the uniform mesh,tie-up of mesh spacing and time interval.It is believed that LBFS has a great potential for solving various flow problems in practice.

    [1] Roach P J.Computational fluid dynamics[M].Hermosa Beach,USA:Hermosa Press,1972.

    [2] Anderson D A,Tannehill J C,Pletcher R H.Computational fluid mechanics and heat transfer[M].New York,USA:McGraw-Hill,1984.

    [3] Hirsch C.Numerical computation of internal and external flows[M].Hoboken,USA:John Wiley &Sons,1988.

    [4] Fletcher C A J.Computational techniques for fluid dynamics:fundamental and general techniques[M].Berlin,Germany:Springer-Verlag,1991.

    [5] Anderson J D.Computational fluid dynamics:the basics with applications[M].New York,USA:McGraw-Hill,1995.

    [6] Versteeg H K,Malalasekera W.An introduction to computational fluid dynamics:the finite volume method[M].Harlow,England:Longman Scientific&Technical,1995.

    [7] Donea J,Huerta A.Finite element methods for flow problems[M].Hoboken,USA:John Wiley,2003.

    [8] Wendt J F.Computational fluid dynamics[M].Berlin,Germany:Springer Berlin Heidelberg,2009.

    [9] Funaro D.Polynomial approximation of differential equations[M].Berlin,Germany:Springer-Verlag,1992.

    [10]Buhmann M D.Radial basis functions:theory and implementations[M].Cambridge University Press,2003.

    [11]Godunov S K.A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[J].Matematicheskii Sbornik,1959,47:271-306.

    [12]Roe P L.Approximate Riemann solvers,parameter vectors,and difference schemes[J].Journal of Computational Physics,1981,43:357-372.

    [13]Steger J,Warming R.Flux vector splitting of the inviscid gas dynamic equations with applications to finite-difference methods[J].Journal of Computational Physics,1981,40:263-293.

    [14]Shu C W,Osher S.Efficient implementation of es-sentially non-oscillatory shock-capturing scheme[J].Journal of Computational Physics,1988,77:439-471.

    [15]Shu C W.High order weighted essentially non-oscillatory schemes for convection dominated problems[J].SIAM Review,2009,51:82-126.

    [16]B van Leer,Lo M.A discontinuous Galerkin method for diffusion based on recovery[J].Journal of Scientific Computation,2011,46:314-328.

    [17]Xu K.A gas-kinetic BGK scheme for the Navier-Stocks equations and its connection with artificial dissipation and Godunov method[J].J Comput Phys,2001,171:289-335.

    [18]Chen S Z,Xu K,Lee C B,et al,A unified gas kinetic scheme with moving mesh and velocity space adaptation[J].Journal of Computational Physics,2012,231:6643-6664.

    [19]Yang L M,Shu C,Wu J,et al.Circular functionbased gas-kinetic scheme for simulation of inviscid compressible flows[J].J Comput Phy,2013,255:540-557.

    [20]Chen S,Chen H,Martínez D,et al.Lattice Boltzmann model for simulation of magnetohydrodynamics[J].Phys Rev Let,1991,67(27):3776-3779.

    [21]Qian Y H,D′Humières D,Lallemand P.Lattice BGK models for Navier-Stokes equation[J].Europhys Lett,1992,17:479-484.

    [22]Chen S,Doolen G.Lattice Boltzmann method for fluid flows[J].Ann Rev Fluid Mech,1998,30:329-64.

    [23]Mei R,Luo L S,Shyy W.An accurate curved boundary treatment in the lattice Boltzmann method[J].J Comput Phys,1999,155:307-330.

    [24]Guo Z L,Shi B C,Wang N C.Lattice BGK model for incompressible Navier-Stokes equation[J].J Comput Phys,2000,165:288-306.

    [25]Shu C,Chew Y T,Niu X D.Least square-based LBM:a meshless approach for simulation of flows with complex geometry[J].Phys Rev E,2001:64,045701.

    [26]Shu C,Niu X D,Chew Y T.Taylor series expansion-and least square-based lattice Boltzmann method:two-dimensional formulation and its applications[J].Phys Rev E,2002,65:036708.

    [27]Succi S,Mesoscopic modeling of slip motion at fluidsolid interfaces with heterogeneous catalysis[J].Phys Rev Lett,2002,89:064502.

    [28]Shan X,Yuan X F,Chen H.Kinetic theory representation of hydrodynamics:A way beyond the Navi-er-Stokes equation[J].J Fluid Mech,2006,550:413-441.

    [29]Guo Z L,Asinari P,Zheng C G.Lattice Boltzmann equation for microscale gas flows of binary mixtures[J].Phys Rev E,2009,79:026702.

    [30]Aidun C K,Clausen J R.Lattice-Boltzmann method for complex flows[J].Ann Rev Fluid Mech,2010,42:439-72.

    [31]Wu J,Shu C.A solution-adaptive lattice Boltzmann method for two-dimensional incompressible viscous flows[J].J Comput Phys,2011,230:2246-2269.

    [32]Guo Z,Shu C.Lattice Boltzmann method and its applications in engineering[J].World Scientific Publishing,2013.

    [33]Yang L M,Shu C,Wu J.Development and comparative studies of three non-free parameter lattice Boltzmann models for simulation of compressible flows[J].Adv Appl Math Mech,2012,4:454-472.

    [34]Yang L M,Shu C,Wu J.A moment conservationbased non-free parameter compressible lattice Boltzmann model and its application for flux evaluation at cell interface[J].Comput Fluids,2013,79:190-199.

    [35]Stolcis L,Johnston L J.Solution of the Euler equations on unstructured grids for two-dimensional compressible flow[J].Aeronautical Journal,1990,94:181-195.

    [36]Woodward P,Colella P.The numerical simulation of two-dimensional fluid flow with strong shocks[J].Journal of Computational Physics,1984,54:115-173.

    [37]Batina J T.Accuracy of an unstructured-grid upwind-Euler algorithm for the ONERA M6wing[J].J Aircraft,1991,28:397-402.

    [38]Schmitt V,Charpin F.Pressure distributions on the ONERA-M6-wing at transonic Mach numbers,experimental data base for computer program assessment[J].Report of the Fluid Dynamics Panel Working Group 04,1979,AGARD AR:138.

    [39]Ghia U,Chia K N,Shin C T.High-resolutions for incompressible flow using the Navier-Stokes equations:a multigrid method[J].J Comput Phys,1982,48:387-411.

    [40]Fuchs L,Tillmark N.Numerical and experimental study of driven flow in a polar cavity[J].Int J Num Methods in Fluids,1985,5:311-329.

    [41]Shu C,Niu X D,Chew Y T.Taylor series expansion-and least square-based lattice Boltzmann method:two-dimensional formulation and its applications[J].Phys Rev E,2002,65:036708.

    [42]Bouard R,Coutanceau M.The early stage of development of the wake behind an impulsively started cylinder for 40<Re<104[J].J Fluid Mech,1980,101:583-607.

    [43]Niu X D,Chew Y T,Shu C.Simulation of flows around an impulsively started circular by Taylor series expansion-and least squares-based lattice Boltzmann method[J].J Comput Phys,2003,188:176-193.

    韩国高清视频一区二区三区| 人体艺术视频欧美日本| 少妇熟女欧美另类| 国产在视频线精品| 久久久久精品性色| 在线观看www视频免费| 香蕉精品网在线| 丝瓜视频免费看黄片| 一区福利在线观看| 国产成人一区二区在线| 99国产综合亚洲精品| 91国产中文字幕| xxx大片免费视频| 观看美女的网站| 久久国内精品自在自线图片| 韩国av在线不卡| 女人高潮潮喷娇喘18禁视频| 婷婷成人精品国产| 免费看av在线观看网站| 亚洲欧美中文字幕日韩二区| 99九九在线精品视频| 91精品伊人久久大香线蕉| 欧美日韩国产mv在线观看视频| 中文字幕av电影在线播放| 久久综合国产亚洲精品| 午夜久久久在线观看| 精品一区二区三区四区五区乱码 | 久久久精品国产亚洲av高清涩受| 熟女少妇亚洲综合色aaa.| 色婷婷av一区二区三区视频| 国产亚洲午夜精品一区二区久久| 热re99久久国产66热| 久久人人爽人人片av| 丁香六月天网| 欧美+日韩+精品| 90打野战视频偷拍视频| 永久免费av网站大全| 日本av免费视频播放| 七月丁香在线播放| 国产精品秋霞免费鲁丝片| 精品午夜福利在线看| 一区二区三区四区激情视频| 在线 av 中文字幕| 久久精品夜色国产| 中文字幕精品免费在线观看视频| 免费观看性生交大片5| 美女国产视频在线观看| 久久免费观看电影| 精品国产乱码久久久久久男人| 国产一区二区三区综合在线观看| 午夜影院在线不卡| 国产一级毛片在线| 热99国产精品久久久久久7| 观看美女的网站| 2018国产大陆天天弄谢| 国产高清国产精品国产三级| 一区在线观看完整版| 满18在线观看网站| 国产又色又爽无遮挡免| 9热在线视频观看99| 在线观看美女被高潮喷水网站| 性少妇av在线| 亚洲人成电影观看| 国产成人精品久久久久久| 久久女婷五月综合色啪小说| 岛国毛片在线播放| 久久毛片免费看一区二区三区| 久久这里有精品视频免费| 老女人水多毛片| 蜜桃国产av成人99| 三级国产精品片| 九色亚洲精品在线播放| 不卡av一区二区三区| 王馨瑶露胸无遮挡在线观看| 国产成人av激情在线播放| 中文精品一卡2卡3卡4更新| 日韩精品有码人妻一区| 国产男人的电影天堂91| 人人妻人人添人人爽欧美一区卜| av卡一久久| av国产久精品久网站免费入址| 丰满少妇做爰视频| 亚洲av免费高清在线观看| 国产老妇伦熟女老妇高清| 99久国产av精品国产电影| 777久久人妻少妇嫩草av网站| 曰老女人黄片| a 毛片基地| 一本大道久久a久久精品| 春色校园在线视频观看| av电影中文网址| 午夜日韩欧美国产| 亚洲色图综合在线观看| 免费av中文字幕在线| 欧美中文综合在线视频| 欧美中文综合在线视频| 国产片内射在线| 国产激情久久老熟女| 男人舔女人的私密视频| 新久久久久国产一级毛片| 自线自在国产av| 乱人伦中国视频| 久久久久久人妻| 精品亚洲乱码少妇综合久久| 国产亚洲av片在线观看秒播厂| 亚洲美女搞黄在线观看| 久久精品国产鲁丝片午夜精品| 女的被弄到高潮叫床怎么办| 欧美日韩精品网址| 五月天丁香电影| 亚洲五月色婷婷综合| 性色av一级| 伊人久久大香线蕉亚洲五| 国产片特级美女逼逼视频| 99久久精品国产国产毛片| 亚洲综合精品二区| 男人爽女人下面视频在线观看| 国产精品三级大全| 国产视频首页在线观看| 黄色怎么调成土黄色| 久久99热这里只频精品6学生| 日本爱情动作片www.在线观看| 久久热在线av| 夫妻午夜视频| 国产亚洲欧美精品永久| 精品一区二区免费观看| 五月伊人婷婷丁香| 婷婷色麻豆天堂久久| 久久99一区二区三区| 免费女性裸体啪啪无遮挡网站| 亚洲精品美女久久久久99蜜臀 | 国产成人精品久久二区二区91 | 高清在线视频一区二区三区| 少妇精品久久久久久久| 亚洲av在线观看美女高潮| 99热全是精品| 欧美精品亚洲一区二区| 亚洲欧美清纯卡通| 国产精品久久久久久av不卡| 亚洲色图综合在线观看| 免费观看性生交大片5| 欧美日韩视频高清一区二区三区二| 多毛熟女@视频| 大片电影免费在线观看免费| 国产精品秋霞免费鲁丝片| 国产精品无大码| 极品人妻少妇av视频| 欧美日韩亚洲国产一区二区在线观看 | 日本wwww免费看| 亚洲综合色惰| 国产97色在线日韩免费| 夫妻午夜视频| 久久午夜综合久久蜜桃| 国产成人av激情在线播放| 看免费av毛片| xxxhd国产人妻xxx| 国精品久久久久久国模美| 性色avwww在线观看| 老司机亚洲免费影院| 日本欧美国产在线视频| 久久久久精品性色| 国产精品熟女久久久久浪| 国产片特级美女逼逼视频| 99热全是精品| 久久精品久久久久久久性| 久久久精品区二区三区| 蜜桃国产av成人99| 99久久人妻综合| 在线观看一区二区三区激情| 国产免费又黄又爽又色| 欧美在线黄色| av视频免费观看在线观看| 亚洲国产日韩一区二区| 永久免费av网站大全| 狠狠婷婷综合久久久久久88av| 叶爱在线成人免费视频播放| 中文天堂在线官网| 国产欧美日韩综合在线一区二区| 91久久精品国产一区二区三区| 人妻人人澡人人爽人人| 丰满饥渴人妻一区二区三| 国产亚洲av片在线观看秒播厂| kizo精华| 美女高潮到喷水免费观看| 天天躁夜夜躁狠狠久久av| 欧美日韩视频精品一区| 日本vs欧美在线观看视频| 日韩中文字幕视频在线看片| 成人国语在线视频| 亚洲av福利一区| 大片电影免费在线观看免费| 国产成人av激情在线播放| 精品国产一区二区三区久久久樱花| 99国产综合亚洲精品| 精品人妻偷拍中文字幕| 最黄视频免费看| 18禁裸乳无遮挡动漫免费视频| 国产激情久久老熟女| 我的亚洲天堂| 久久午夜福利片| 一区二区三区激情视频| 高清在线视频一区二区三区| 另类亚洲欧美激情| 美女高潮到喷水免费观看| 人体艺术视频欧美日本| 观看av在线不卡| 777米奇影视久久| 久久97久久精品| 国产精品嫩草影院av在线观看| 国产精品嫩草影院av在线观看| 国产精品久久久久久久久免| www.自偷自拍.com| 精品人妻一区二区三区麻豆| 久久婷婷青草| 国产av码专区亚洲av| 一级,二级,三级黄色视频| 日本av免费视频播放| 国产一区亚洲一区在线观看| 国产人伦9x9x在线观看 | 秋霞伦理黄片| 精品人妻熟女毛片av久久网站| 叶爱在线成人免费视频播放| 亚洲精品久久久久久婷婷小说| 国产精品亚洲av一区麻豆 | 国产视频首页在线观看| 国产1区2区3区精品| 天堂俺去俺来也www色官网| 日本91视频免费播放| 另类亚洲欧美激情| 亚洲av电影在线进入| 国产成人精品久久二区二区91 | 巨乳人妻的诱惑在线观看| a级片在线免费高清观看视频| 波野结衣二区三区在线| 1024视频免费在线观看| 搡女人真爽免费视频火全软件| 母亲3免费完整高清在线观看 | 人妻一区二区av| 久久久久国产网址| 精品一区二区三卡| 99国产精品免费福利视频| 最新的欧美精品一区二区| 爱豆传媒免费全集在线观看| 国产av国产精品国产| 人妻系列 视频| 亚洲精华国产精华液的使用体验| 日本-黄色视频高清免费观看| 国产黄色免费在线视频| 18禁裸乳无遮挡动漫免费视频| 校园人妻丝袜中文字幕| 亚洲精品日本国产第一区| 国产欧美日韩一区二区三区在线| 久久这里有精品视频免费| 国产男女超爽视频在线观看| 大香蕉久久成人网| 日韩大片免费观看网站| 久久综合国产亚洲精品| 国产精品久久久久久精品电影小说| 另类精品久久| 国产xxxxx性猛交| 国产免费现黄频在线看| 欧美日韩一级在线毛片| 日韩成人av中文字幕在线观看| a级毛片黄视频| 丝瓜视频免费看黄片| 秋霞伦理黄片| 日韩av在线免费看完整版不卡| 视频区图区小说| 热re99久久国产66热| 国产一区有黄有色的免费视频| 亚洲一区二区三区欧美精品| 啦啦啦在线免费观看视频4| 亚洲图色成人| 永久网站在线| 黄片播放在线免费| 欧美亚洲日本最大视频资源| 国产伦理片在线播放av一区| 免费高清在线观看视频在线观看| 91久久精品国产一区二区三区| 在线 av 中文字幕| 久久久久久人妻| 18禁国产床啪视频网站| 成人毛片60女人毛片免费| a级片在线免费高清观看视频| 亚洲精品美女久久av网站| 91午夜精品亚洲一区二区三区| 日韩不卡一区二区三区视频在线| 亚洲av电影在线观看一区二区三区| 国产在视频线精品| 欧美成人精品欧美一级黄| 久热这里只有精品99| 国产免费现黄频在线看| 国产极品天堂在线| av片东京热男人的天堂| 久久影院123| 人妻人人澡人人爽人人| 妹子高潮喷水视频| 午夜福利一区二区在线看| 欧美激情高清一区二区三区 | 日韩av在线免费看完整版不卡| 各种免费的搞黄视频| 久久久久久久精品精品| 国产老妇伦熟女老妇高清| 欧美成人午夜免费资源| 欧美激情 高清一区二区三区| 亚洲一码二码三码区别大吗| 欧美日韩一级在线毛片| 久久ye,这里只有精品| 一区二区三区乱码不卡18| 麻豆精品久久久久久蜜桃| 国产精品久久久av美女十八| 岛国毛片在线播放| 亚洲色图 男人天堂 中文字幕| 日韩三级伦理在线观看| 啦啦啦在线免费观看视频4| 久久精品夜色国产| 欧美在线黄色| 亚洲中文av在线| 欧美日韩成人在线一区二区| 日韩制服骚丝袜av| 欧美精品国产亚洲| 毛片一级片免费看久久久久| 国产精品久久久久久精品古装| 午夜福利一区二区在线看| 国产福利在线免费观看视频| 少妇的丰满在线观看| 卡戴珊不雅视频在线播放| 国产精品99久久99久久久不卡 | 91精品伊人久久大香线蕉| 国产精品一二三区在线看| 国产xxxxx性猛交| 午夜福利在线免费观看网站| 一区二区av电影网| 亚洲精品av麻豆狂野| 999精品在线视频| 亚洲第一青青草原| 国产黄色免费在线视频| 丰满乱子伦码专区| a 毛片基地| 国产老妇伦熟女老妇高清| 午夜影院在线不卡| 中文字幕人妻丝袜制服| 性少妇av在线| 亚洲精品国产av成人精品| 欧美av亚洲av综合av国产av | 伦精品一区二区三区| 免费人妻精品一区二区三区视频| 国产日韩一区二区三区精品不卡| 侵犯人妻中文字幕一二三四区| 成年动漫av网址| 伊人亚洲综合成人网| 久久热在线av| 如何舔出高潮| 成人亚洲精品一区在线观看| 免费观看a级毛片全部| 卡戴珊不雅视频在线播放| 国产亚洲一区二区精品| 久久精品国产亚洲av涩爱| 99热国产这里只有精品6| 日韩欧美精品免费久久| 久久久a久久爽久久v久久| 午夜福利网站1000一区二区三区| 久久99热这里只频精品6学生| 免费黄频网站在线观看国产| 亚洲第一区二区三区不卡| 亚洲国产欧美在线一区| 久久久久精品性色| 久久久久精品人妻al黑| 七月丁香在线播放| 亚洲三区欧美一区| 国产白丝娇喘喷水9色精品| 一级毛片电影观看| 97精品久久久久久久久久精品| 欧美激情 高清一区二区三区| 日本av免费视频播放| 91在线精品国自产拍蜜月| 99久久人妻综合| 叶爱在线成人免费视频播放| 中文精品一卡2卡3卡4更新| 久久久久国产精品人妻一区二区| 国产毛片在线视频| 中文字幕人妻熟女乱码| 国产精品久久久久久av不卡| av在线老鸭窝| 天天躁夜夜躁狠狠久久av| 搡老乐熟女国产| 狠狠精品人妻久久久久久综合| 欧美最新免费一区二区三区| 亚洲综合色惰| 国产精品久久久久成人av| 天天影视国产精品| 男女边摸边吃奶| 亚洲情色 制服丝袜| 国产欧美日韩一区二区三区在线| videosex国产| 各种免费的搞黄视频| 男人爽女人下面视频在线观看| 国产精品国产av在线观看| 国产精品无大码| 嫩草影院入口| 亚洲国产毛片av蜜桃av| 伊人亚洲综合成人网| 久久影院123| 亚洲国产成人一精品久久久| 国产白丝娇喘喷水9色精品| av片东京热男人的天堂| 欧美日韩av久久| 久久av网站| 看非洲黑人一级黄片| videossex国产| 亚洲三级黄色毛片| av在线观看视频网站免费| 免费观看无遮挡的男女| videos熟女内射| 久久99一区二区三区| 久久久久人妻精品一区果冻| 丝袜在线中文字幕| 大陆偷拍与自拍| 午夜激情久久久久久久| 亚洲欧美日韩另类电影网站| 九草在线视频观看| 18禁裸乳无遮挡动漫免费视频| 9热在线视频观看99| 欧美+日韩+精品| 女人精品久久久久毛片| 国产精品99久久99久久久不卡 | av在线老鸭窝| 亚洲av国产av综合av卡| 我的亚洲天堂| 巨乳人妻的诱惑在线观看| 91国产中文字幕| 亚洲婷婷狠狠爱综合网| 精品视频人人做人人爽| 两性夫妻黄色片| 一级毛片我不卡| 人妻 亚洲 视频| 日本午夜av视频| 日本欧美国产在线视频| 国产日韩欧美在线精品| 国产福利在线免费观看视频| 午夜免费鲁丝| 精品亚洲成a人片在线观看| 日韩人妻精品一区2区三区| 黄色怎么调成土黄色| 国产色婷婷99| 亚洲熟女精品中文字幕| 亚洲四区av| 午夜福利,免费看| 一级毛片 在线播放| 国产乱来视频区| 丝袜美腿诱惑在线| 伦理电影免费视频| 伊人久久大香线蕉亚洲五| 亚洲av日韩在线播放| 国产精品女同一区二区软件| 在线观看人妻少妇| 成人黄色视频免费在线看| 一本大道久久a久久精品| 国产成人一区二区在线| 中文字幕色久视频| 婷婷色av中文字幕| 熟女电影av网| 日韩成人av中文字幕在线观看| 欧美av亚洲av综合av国产av | 纵有疾风起免费观看全集完整版| 精品少妇久久久久久888优播| av不卡在线播放| 天天躁狠狠躁夜夜躁狠狠躁| 亚洲激情五月婷婷啪啪| 少妇猛男粗大的猛烈进出视频| 交换朋友夫妻互换小说| 亚洲欧美日韩另类电影网站| 91精品三级在线观看| 十分钟在线观看高清视频www| 韩国高清视频一区二区三区| 精品国产乱码久久久久久男人| 91精品国产国语对白视频| 丝瓜视频免费看黄片| 国产黄色免费在线视频| 亚洲成人手机| 18禁国产床啪视频网站| 一级毛片 在线播放| 欧美亚洲日本最大视频资源| 亚洲欧美精品自产自拍| av国产精品久久久久影院| 男女无遮挡免费网站观看| 日本午夜av视频| 水蜜桃什么品种好| 黑丝袜美女国产一区| 两个人看的免费小视频| 国产又爽黄色视频| 最黄视频免费看| 免费黄色在线免费观看| 免费观看av网站的网址| 久久97久久精品| 国产又色又爽无遮挡免| 捣出白浆h1v1| 午夜av观看不卡| 少妇 在线观看| 日本-黄色视频高清免费观看| 最近中文字幕高清免费大全6| 国产色婷婷99| 精品一区二区免费观看| 成人二区视频| 亚洲精品日韩在线中文字幕| 美女高潮到喷水免费观看| 多毛熟女@视频| 亚洲国产毛片av蜜桃av| 最新中文字幕久久久久| 亚洲三区欧美一区| 久久精品国产亚洲av高清一级| 91精品国产国语对白视频| 亚洲精华国产精华液的使用体验| 一本色道久久久久久精品综合| 午夜精品国产一区二区电影| 国产精品国产三级专区第一集| 日本欧美国产在线视频| 女性生殖器流出的白浆| 久久久久人妻精品一区果冻| 日韩视频在线欧美| 亚洲av男天堂| 亚洲精品国产av蜜桃| 久久99精品国语久久久| 国产成人aa在线观看| 久久精品国产亚洲av高清一级| 亚洲婷婷狠狠爱综合网| 看免费成人av毛片| 久久午夜综合久久蜜桃| 久久婷婷青草| 各种免费的搞黄视频| 久久精品国产a三级三级三级| 亚洲天堂av无毛| 久久久久久免费高清国产稀缺| 2022亚洲国产成人精品| 一二三四在线观看免费中文在| 久久久a久久爽久久v久久| 又大又黄又爽视频免费| 丝袜在线中文字幕| 黄片无遮挡物在线观看| 国产乱来视频区| 视频在线观看一区二区三区| 欧美日韩av久久| 永久网站在线| 99久国产av精品国产电影| 精品国产露脸久久av麻豆| 午夜福利视频精品| 黄色一级大片看看| 午夜福利在线观看免费完整高清在| 色视频在线一区二区三区| 熟女av电影| 久久久久人妻精品一区果冻| 少妇的逼水好多| 国产片内射在线| 80岁老熟妇乱子伦牲交| 多毛熟女@视频| 人妻 亚洲 视频| 黑人猛操日本美女一级片| 一区二区日韩欧美中文字幕| 有码 亚洲区| 国产精品香港三级国产av潘金莲 | 亚洲精品国产av蜜桃| 26uuu在线亚洲综合色| 亚洲成人一二三区av| av在线观看视频网站免费| 侵犯人妻中文字幕一二三四区| 精品人妻偷拍中文字幕| 亚洲国产精品成人久久小说| 久久精品国产亚洲av高清一级| 中文字幕精品免费在线观看视频| 免费高清在线观看日韩| 国产免费一区二区三区四区乱码| videosex国产| 精品少妇内射三级| 性高湖久久久久久久久免费观看| av国产精品久久久久影院| 在线天堂中文资源库| 校园人妻丝袜中文字幕| 99re6热这里在线精品视频| 黑人欧美特级aaaaaa片| 午夜福利,免费看| 黄片播放在线免费| 高清视频免费观看一区二区| 久久久久久久久免费视频了| 欧美黄色片欧美黄色片| 国产精品女同一区二区软件| 亚洲,欧美,日韩| 日本爱情动作片www.在线观看| 免费在线观看完整版高清| 黄片播放在线免费| 深夜精品福利| 日韩成人av中文字幕在线观看| 久久精品国产亚洲av涩爱| 亚洲精品久久午夜乱码| 日韩,欧美,国产一区二区三区| 久久青草综合色| 国产在线免费精品| 毛片一级片免费看久久久久| 免费高清在线观看日韩| 人妻人人澡人人爽人人| 狠狠婷婷综合久久久久久88av| av不卡在线播放| 男女下面插进去视频免费观看| 欧美日韩精品网址| 这个男人来自地球电影免费观看 | av不卡在线播放| 高清黄色对白视频在线免费看| a级毛片黄视频| 91午夜精品亚洲一区二区三区| 免费黄色在线免费观看| 欧美日韩成人在线一区二区| 国产精品久久久久久精品古装| 制服诱惑二区| 免费黄网站久久成人精品| 亚洲伊人久久精品综合| tube8黄色片|