• <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.www免费av| 亚洲国产欧美网| 免费av不卡在线播放| 国产精品香港三级国产av潘金莲| 成在线人永久免费视频| 亚洲国产色片| 99精品在免费线老司机午夜| 国产极品精品免费视频能看的| 一区福利在线观看| 国产精品 国内视频| 久久精品国产99精品国产亚洲性色| 日日摸夜夜添夜夜添小说| 夜夜看夜夜爽夜夜摸| 国内揄拍国产精品人妻在线| 国产成人精品久久二区二区91| 女人被狂操c到高潮| 看黄色毛片网站| 亚洲国产高清在线一区二区三| 久久精品aⅴ一区二区三区四区| 亚洲欧美一区二区三区黑人| 亚洲男人的天堂狠狠| 韩国av一区二区三区四区| 亚洲 欧美 日韩 在线 免费| 深夜精品福利| 俺也久久电影网| 国产精品久久电影中文字幕| av福利片在线观看| 国产精品久久久久久精品电影| cao死你这个sao货| 欧美日韩亚洲国产一区二区在线观看| 国产乱人伦免费视频| 五月伊人婷婷丁香| 欧美黑人巨大hd| 欧美另类亚洲清纯唯美| av天堂中文字幕网| 1024香蕉在线观看| 国产精品 国内视频| 免费搜索国产男女视频| 亚洲欧美精品综合一区二区三区| 欧美乱色亚洲激情| 国产毛片a区久久久久| 中文亚洲av片在线观看爽| 欧美三级亚洲精品| 国产av在哪里看| 亚洲片人在线观看| 亚洲精品美女久久av网站| 亚洲熟妇熟女久久| 欧美成人性av电影在线观看| 亚洲五月天丁香| 在线观看午夜福利视频| 国产真实乱freesex| 在线免费观看不下载黄p国产 | 日本在线视频免费播放| 怎么达到女性高潮| 老司机深夜福利视频在线观看| 亚洲人成网站在线播放欧美日韩| 嫩草影视91久久| 亚洲片人在线观看| 久久久久久国产a免费观看| 69av精品久久久久久| 国产真实乱freesex| 欧美中文综合在线视频| 丁香六月欧美| 大型黄色视频在线免费观看| 偷拍熟女少妇极品色| 亚洲人成伊人成综合网2020| 国产免费男女视频| 国产伦精品一区二区三区视频9 | 黄片大片在线免费观看| 观看免费一级毛片| 久久精品91蜜桃| 亚洲精品色激情综合| 国产精华一区二区三区| 小蜜桃在线观看免费完整版高清| 少妇人妻一区二区三区视频| 成人无遮挡网站| 亚洲天堂国产精品一区在线| ponron亚洲| 久久久久久国产a免费观看| 热99re8久久精品国产| 在线国产一区二区在线| 日韩免费av在线播放| 久久久久久久午夜电影| 久久久久国产一级毛片高清牌| 国产毛片a区久久久久| 亚洲美女视频黄频| 久久午夜亚洲精品久久| 亚洲av片天天在线观看| 中文字幕人成人乱码亚洲影| 久久天躁狠狠躁夜夜2o2o| 黄色丝袜av网址大全| 在线a可以看的网站| 色吧在线观看| 男人和女人高潮做爰伦理| 亚洲 欧美一区二区三区| 又粗又爽又猛毛片免费看| 综合色av麻豆| 亚洲精品456在线播放app | 欧美成狂野欧美在线观看| 91老司机精品| 久久久久久九九精品二区国产| 久久久久久国产a免费观看| 在线看三级毛片| 熟妇人妻久久中文字幕3abv| 国产91精品成人一区二区三区| 两性夫妻黄色片| 99精品在免费线老司机午夜| 国产一区二区在线观看日韩 | 操出白浆在线播放| 老司机深夜福利视频在线观看| 国内揄拍国产精品人妻在线| 最近在线观看免费完整版| 波多野结衣巨乳人妻| 亚洲乱码一区二区免费版| 日日干狠狠操夜夜爽| www.熟女人妻精品国产| 91字幕亚洲| 国产精品女同一区二区软件 | 99精品欧美一区二区三区四区| 国产淫片久久久久久久久 | 亚洲欧美激情综合另类| 亚洲va日本ⅴa欧美va伊人久久| 热99在线观看视频| 成人无遮挡网站| 久久热在线av| 身体一侧抽搐| 欧美日韩中文字幕国产精品一区二区三区| 97超级碰碰碰精品色视频在线观看| 国产免费男女视频| 亚洲av免费在线观看| 99久久国产精品久久久| 国产精品一区二区三区四区免费观看 | av国产免费在线观看| 亚洲国产欧美一区二区综合| 小蜜桃在线观看免费完整版高清| 国产真人三级小视频在线观看| 又大又爽又粗| netflix在线观看网站| 精华霜和精华液先用哪个| 久久精品91无色码中文字幕| 高潮久久久久久久久久久不卡| 欧美+亚洲+日韩+国产| 国产精品女同一区二区软件 | 亚洲国产欧美人成| 精品久久久久久成人av| 日韩人妻高清精品专区| 一个人观看的视频www高清免费观看 | 黄色 视频免费看| or卡值多少钱| 日韩欧美在线二视频| 欧美三级亚洲精品| 搡老妇女老女人老熟妇| 一级毛片女人18水好多| 91av网一区二区| 18美女黄网站色大片免费观看| 国内精品久久久久精免费| 亚洲精品在线观看二区| 18禁裸乳无遮挡免费网站照片| 午夜精品一区二区三区免费看| 精品人妻1区二区| 色播亚洲综合网| 天堂网av新在线| 亚洲精品一区av在线观看| 免费人成视频x8x8入口观看| 99国产精品99久久久久| 国产免费av片在线观看野外av| 后天国语完整版免费观看| 黑人操中国人逼视频| 午夜福利欧美成人| 国产麻豆成人av免费视频| 亚洲最大成人中文| 老汉色av国产亚洲站长工具| 精品日产1卡2卡| 亚洲 欧美 日韩 在线 免费| 日日干狠狠操夜夜爽| 国产精品一区二区免费欧美| 五月伊人婷婷丁香| 狂野欧美白嫩少妇大欣赏| 欧美日韩瑟瑟在线播放| 欧美一区二区国产精品久久精品| 欧美中文综合在线视频| 免费看a级黄色片| 嫩草影院入口| 香蕉久久夜色| 亚洲熟妇中文字幕五十中出| 亚洲国产精品久久男人天堂| 亚洲精品一区av在线观看| 久久精品91无色码中文字幕| 熟妇人妻久久中文字幕3abv| 亚洲精品国产精品久久久不卡| 91麻豆av在线| 大型黄色视频在线免费观看| 日本一二三区视频观看| 成人高潮视频无遮挡免费网站| 99热只有精品国产| 午夜福利成人在线免费观看| 午夜福利在线在线| 免费观看的影片在线观看| 久久精品国产综合久久久| 九九久久精品国产亚洲av麻豆 | 九九在线视频观看精品| 黄色视频,在线免费观看| 亚洲第一电影网av| 精品久久久久久,| 99riav亚洲国产免费| 日韩欧美在线乱码| 精品国产美女av久久久久小说| 黄色片一级片一级黄色片| 国语自产精品视频在线第100页| 九色成人免费人妻av| 99久国产av精品| 身体一侧抽搐| 欧美成人一区二区免费高清观看 | 久久久久久大精品| 一进一出抽搐动态| 特级一级黄色大片| 国产精品久久久久久精品电影| 中国美女看黄片| 欧美大码av| 黄色日韩在线| 男女下面进入的视频免费午夜| 91在线精品国自产拍蜜月 | svipshipincom国产片| 国产乱人伦免费视频| 午夜免费观看网址| 久久精品91无色码中文字幕| 亚洲第一电影网av| 国产一区二区三区视频了| 中文字幕人妻丝袜一区二区| 亚洲色图 男人天堂 中文字幕| 99热精品在线国产| 日韩人妻高清精品专区| 极品教师在线免费播放| 亚洲一区高清亚洲精品| www日本在线高清视频| 久久精品aⅴ一区二区三区四区| 亚洲欧美一区二区三区黑人| 亚洲欧洲精品一区二区精品久久久| 亚洲 国产 在线| 一边摸一边抽搐一进一小说| 岛国在线观看网站| 久久中文字幕一级| 色综合婷婷激情| 18禁国产床啪视频网站| 亚洲avbb在线观看| 村上凉子中文字幕在线| 三级毛片av免费| 熟女电影av网| 久久天躁狠狠躁夜夜2o2o| 搡老熟女国产l中国老女人| 欧美黑人欧美精品刺激| 99久久成人亚洲精品观看| 日本与韩国留学比较| 欧美乱妇无乱码| 国产精品乱码一区二三区的特点| 国产精品日韩av在线免费观看| 级片在线观看| 美女免费视频网站| av女优亚洲男人天堂 | 给我免费播放毛片高清在线观看| 三级男女做爰猛烈吃奶摸视频| 亚洲九九香蕉| 亚洲av成人一区二区三| 深夜精品福利| 国产 一区 欧美 日韩| 国产成人系列免费观看| 97超级碰碰碰精品色视频在线观看| 午夜福利在线观看免费完整高清在 | 禁无遮挡网站| 国产精品1区2区在线观看.| 舔av片在线| 午夜精品一区二区三区免费看| 久久亚洲精品不卡| 特大巨黑吊av在线直播| 久久国产乱子伦精品免费另类| 成人性生交大片免费视频hd| av在线蜜桃| 亚洲精华国产精华精| 亚洲自偷自拍图片 自拍| 国产精品九九99| 99久久成人亚洲精品观看| 成人三级黄色视频| 国产亚洲欧美98| 国内少妇人妻偷人精品xxx网站 | 欧美av亚洲av综合av国产av| 国产亚洲欧美98| 老熟妇仑乱视频hdxx| 国产精品久久视频播放| 婷婷精品国产亚洲av| 久久中文字幕一级| 亚洲国产日韩欧美精品在线观看 | 亚洲人与动物交配视频| 三级国产精品欧美在线观看 | 欧美一区二区精品小视频在线| 18美女黄网站色大片免费观看| 香蕉久久夜色| 久久中文字幕一级| 国产午夜精品论理片| 麻豆国产97在线/欧美| 99久久成人亚洲精品观看| 日韩 欧美 亚洲 中文字幕| av片东京热男人的天堂| 超碰成人久久| 美女黄网站色视频| 国产亚洲欧美98| 欧美日韩综合久久久久久 | 听说在线观看完整版免费高清| 国产人伦9x9x在线观看| 欧美又色又爽又黄视频| 国产精品久久视频播放| 制服人妻中文乱码| 1000部很黄的大片| 久久久精品欧美日韩精品| 91麻豆av在线| 动漫黄色视频在线观看| 国产精品爽爽va在线观看网站| 精品日产1卡2卡| 日本黄色视频三级网站网址| 女生性感内裤真人,穿戴方法视频| 国产单亲对白刺激| 亚洲专区国产一区二区| 九色成人免费人妻av| 中出人妻视频一区二区| 国产精品久久久久久久电影 | 国产成人精品久久二区二区91| xxx96com| av中文乱码字幕在线| 久久天躁狠狠躁夜夜2o2o| 成人永久免费在线观看视频| 国产精品女同一区二区软件 | 国产精品久久视频播放| 又紧又爽又黄一区二区| 一本久久中文字幕| 成年女人看的毛片在线观看| 男人和女人高潮做爰伦理| 99视频精品全部免费 在线 | 精品人妻1区二区| 国产成人av教育| 国产爱豆传媒在线观看| 国产成人一区二区三区免费视频网站| ponron亚洲| av女优亚洲男人天堂 | 首页视频小说图片口味搜索| 久久国产乱子伦精品免费另类| 欧美色欧美亚洲另类二区| 国产乱人视频| 精品一区二区三区四区五区乱码| 神马国产精品三级电影在线观看| 免费在线观看亚洲国产| 久久久久亚洲av毛片大全| 女警被强在线播放| 国产精品影院久久| e午夜精品久久久久久久| 国产av在哪里看| 老司机午夜十八禁免费视频| 2021天堂中文幕一二区在线观| 亚洲欧美激情综合另类| 亚洲自拍偷在线| 日本撒尿小便嘘嘘汇集6| 婷婷精品国产亚洲av在线| 一进一出抽搐动态| 日本五十路高清| 午夜激情欧美在线| 亚洲乱码一区二区免费版| 夜夜夜夜夜久久久久| 亚洲精品在线观看二区| 欧美中文日本在线观看视频| 男人舔女人下体高潮全视频| 蜜桃久久精品国产亚洲av| 最新中文字幕久久久久 | 午夜免费成人在线视频| 女人被狂操c到高潮| 亚洲欧美日韩无卡精品| 成熟少妇高潮喷水视频| 日韩欧美一区二区三区在线观看| av黄色大香蕉| 国产成年人精品一区二区| 国产成人av激情在线播放| 波多野结衣高清作品| 久久人妻av系列| 宅男免费午夜| 亚洲欧美日韩高清在线视频| 麻豆国产97在线/欧美| 国产av麻豆久久久久久久| 久久香蕉国产精品| 在线观看舔阴道视频| 国产精品av久久久久免费| 国产亚洲精品一区二区www| 俄罗斯特黄特色一大片| 又黄又爽又免费观看的视频| 免费看十八禁软件| 国产极品精品免费视频能看的| 久久久色成人| 欧美+亚洲+日韩+国产| 精品一区二区三区四区五区乱码| 少妇人妻一区二区三区视频| 法律面前人人平等表现在哪些方面| 黄色女人牲交| 夜夜躁狠狠躁天天躁| 99在线人妻在线中文字幕| 亚洲av中文字字幕乱码综合| 偷拍熟女少妇极品色| 日韩国内少妇激情av| 69av精品久久久久久| 国产精品1区2区在线观看.| 这个男人来自地球电影免费观看| 亚洲一区二区三区色噜噜| 亚洲欧美日韩高清专用| 此物有八面人人有两片| 又爽又黄无遮挡网站| 欧美在线黄色| 国产久久久一区二区三区| 国产精品久久久人人做人人爽| 最新美女视频免费是黄的| 日本精品一区二区三区蜜桃| 麻豆一二三区av精品| 国产精品爽爽va在线观看网站| 欧美日韩国产亚洲二区| 狠狠狠狠99中文字幕| 激情在线观看视频在线高清| 高清毛片免费观看视频网站| 悠悠久久av| 非洲黑人性xxxx精品又粗又长| 欧美日本亚洲视频在线播放| 欧美成人免费av一区二区三区| 热99在线观看视频| 亚洲人成网站高清观看| 亚洲无线在线观看| 在线十欧美十亚洲十日本专区| 久久亚洲真实| 身体一侧抽搐| 97超级碰碰碰精品色视频在线观看| 搡老岳熟女国产| 亚洲男人的天堂狠狠| 国产私拍福利视频在线观看| 窝窝影院91人妻| 黄色视频,在线免费观看| 757午夜福利合集在线观看| 亚洲精华国产精华精| 成年女人看的毛片在线观看| 99热6这里只有精品| 国产成人系列免费观看| 天堂影院成人在线观看| 成年版毛片免费区| 免费搜索国产男女视频| 亚洲中文字幕一区二区三区有码在线看 | 国产v大片淫在线免费观看| 亚洲精品456在线播放app | 久久精品亚洲精品国产色婷小说| 国产精品一及| 亚洲成a人片在线一区二区| 亚洲成人久久性| 精品久久久久久久末码| 国产午夜精品论理片| 国产午夜精品久久久久久| 亚洲天堂国产精品一区在线| 欧美在线黄色| 欧美黄色淫秽网站| 999久久久国产精品视频| 国产亚洲精品久久久com| 免费无遮挡裸体视频| 国内毛片毛片毛片毛片毛片| 免费在线观看日本一区| 九九久久精品国产亚洲av麻豆 | 搞女人的毛片| 淫秽高清视频在线观看| 91老司机精品| 国产成年人精品一区二区| 久久人妻av系列| 亚洲中文av在线| 不卡av一区二区三区| 又大又爽又粗| 欧美又色又爽又黄视频| 午夜成年电影在线免费观看| 国产午夜精品久久久久久| 一个人看视频在线观看www免费 | 亚洲第一电影网av| 亚洲国产欧美人成| 日韩精品中文字幕看吧| 给我免费播放毛片高清在线观看| 国产av在哪里看| 日本三级黄在线观看| 99在线人妻在线中文字幕| 久久精品亚洲精品国产色婷小说| 国产成人影院久久av| 欧美+亚洲+日韩+国产| 成人亚洲精品av一区二区| 国产成年人精品一区二区| 亚洲成av人片在线播放无| 99久久无色码亚洲精品果冻| 三级国产精品欧美在线观看 | 宅男免费午夜| 亚洲成av人片在线播放无| 久久亚洲真实| 日韩免费av在线播放| 色综合欧美亚洲国产小说| 欧美黄色片欧美黄色片| 欧美日韩国产亚洲二区| 高清毛片免费观看视频网站| 国产av在哪里看| 丰满人妻熟妇乱又伦精品不卡| 亚洲精品国产精品久久久不卡| 成人永久免费在线观看视频| 黄片大片在线免费观看| 国产三级在线视频| 国产 一区 欧美 日韩| 嫩草影院精品99| 亚洲专区国产一区二区| 精品无人区乱码1区二区| 亚洲美女视频黄频| 亚洲av中文字字幕乱码综合| 男女视频在线观看网站免费| 男人舔女人的私密视频| 欧洲精品卡2卡3卡4卡5卡区| 久久精品综合一区二区三区| 一个人看的www免费观看视频| 成人性生交大片免费视频hd| 少妇丰满av| 欧美最黄视频在线播放免费| 国产成人精品无人区| 麻豆久久精品国产亚洲av| 老司机午夜十八禁免费视频| 亚洲激情在线av| 久久久久久久久中文| a级毛片在线看网站| 好男人电影高清在线观看| 法律面前人人平等表现在哪些方面| 亚洲精品粉嫩美女一区| 亚洲午夜精品一区,二区,三区| 一级作爱视频免费观看| 男女下面进入的视频免费午夜| 成人一区二区视频在线观看| av视频在线观看入口| 成人av一区二区三区在线看| 成熟少妇高潮喷水视频| 国产一区二区在线av高清观看| 欧美不卡视频在线免费观看| 一进一出抽搐动态| 国产精品久久久人人做人人爽| 在线国产一区二区在线| 十八禁人妻一区二区| 村上凉子中文字幕在线| 久久久精品大字幕| 黄频高清免费视频| 每晚都被弄得嗷嗷叫到高潮| 成人亚洲精品av一区二区| 床上黄色一级片| 免费无遮挡裸体视频| 国产麻豆成人av免费视频| 男女做爰动态图高潮gif福利片| 真人做人爱边吃奶动态| 免费看十八禁软件| 午夜福利18| 夜夜躁狠狠躁天天躁| 欧美成人性av电影在线观看| 999久久久国产精品视频| 欧美国产日韩亚洲一区| 亚洲性夜色夜夜综合| 亚洲国产高清在线一区二区三| 欧美中文日本在线观看视频| 亚洲精品美女久久久久99蜜臀| 亚洲国产精品合色在线| 日韩有码中文字幕| 男女做爰动态图高潮gif福利片| 欧美高清成人免费视频www| 俺也久久电影网| 日韩欧美精品v在线| 精品久久久久久久久久免费视频| 两个人视频免费观看高清| 悠悠久久av| 日韩欧美在线二视频| 最近视频中文字幕2019在线8| 一区二区三区激情视频| 亚洲欧美精品综合久久99| 亚洲国产欧美人成| 黑人巨大精品欧美一区二区mp4| 欧美丝袜亚洲另类 | 韩国av一区二区三区四区| 亚洲成av人片免费观看| 国产又色又爽无遮挡免费看| 人人妻人人看人人澡| 桃红色精品国产亚洲av| 一级a爱片免费观看的视频| 欧美最黄视频在线播放免费| 色精品久久人妻99蜜桃| 亚洲自拍偷在线| 窝窝影院91人妻| 一边摸一边抽搐一进一小说| 美女大奶头视频| 日日摸夜夜添夜夜添小说| 这个男人来自地球电影免费观看| 少妇人妻一区二区三区视频| 99久久精品一区二区三区| 夜夜爽天天搞| 国产野战对白在线观看| 成人一区二区视频在线观看| 国模一区二区三区四区视频 | 久久人妻av系列| 亚洲av熟女| 波多野结衣高清无吗| 精品一区二区三区av网在线观看| 国产精品久久久久久人妻精品电影| 久久欧美精品欧美久久欧美| 91老司机精品| 国产精品 国内视频| 亚洲自拍偷在线| 国产激情偷乱视频一区二区| a在线观看视频网站| 美女午夜性视频免费| 午夜日韩欧美国产| 久久久色成人| 久久久成人免费电影| 欧美乱色亚洲激情| 18禁黄网站禁片午夜丰满| 久久久久久久久久黄片|