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

    Simulations of Blood Drop Spreading and Impact for Bloodstain Pattern Analysis

    2014-04-17 11:06:05ChuWangandLucyZhang

    Chu Wang and Lucy T.Zhang

    1 Introduction

    The bloodstain pattern analysis(BPA),which is a part of the forensic science,can be crucial in solving a crime scene[Karger,Rand,Fracasso,and Pfeiffer(2008)].It provides more comprehensive evidence for criminologists and scientists by an-alyzing blood generation,motion and spatter patterns.In the most recent review article on fluid dynamics and bloodstain pattern analysis,Attinger et al.[Attinger,Moore,Donaldson,Jafari,and Stone(2013)]addressed the essential role fluid dynamics plays in forensic science.In both communities,BPA and fluid dynamics,the blood drop dynamics as it interacts with surrounding environments,air and/or other liquids or solids,has always been one of the key components of the analysis.The underlying physics of these multi-phase interactions involve blood drop trajectories in air(multi- fluid flows)and the spatter patterns(blood spatter analysis)as blood drops impact on solid surfaces( fluid-solid interactions).Therefore, fluid dynamic analysis,focusing on multiphase flows such as the studies of drop impact and spreading,can be used to perform BPA considering the similarities drawn between the two[Attinger,Moore,Donaldson,Jafari,and Stone(2013)].

    Among the many possible dominant factors such as drop generation,trajectories,impact,staining patterns on different types of surfaces in analyzing blood drops and blood stains,here we will only focus our effort on two of them:(1)drop spreading,and(2)drop impact.Numerical model and simulations of blood drop spreading and impact can offer scientists a visual and dynamic prediction on the exact sequence of events.However,the complicated behaviors of blood drops require the numerical models to be capable of accurately capturing and predicting the intricate interplay among the interfacial forces,which include:(a)the three-phase(air-liquid-solid)dynamic contact line,(b)solving the Navier-Stokes equations with properly modeled surface tension force,and(c)coping with geometry topology changes of the blood drops during the process,such as blood spattering due to impact.

    In the forefront of modelingdrop spreading,which involves the overall surface energy balance among air-liquid,liquid-surface,and air-surface,several numerical algorithms have been proposed and widely used,such as the front tracking method[Tryggvason,Bunner,Esmaeeli,Juric,Al-Rawahi,Tauber,Han,Nas,and Jan(2001);Unverdi and Tryggvason(1992);Witteveen,Koren,and Bakker(2007);Muradoglu and Tryggvason(2008)],volume of fluid method(VOF)[Hirt and Nichols(1981);Cervone,Manservisi,and Scardovelli(2010)]and the level-set method[Xu,Li,Lowengrub,and Zhao(2006);Herrmann(2008)].The classic usage of these methods is only to model air-liquid multiphase flows,which sets the limitations of their applicability.Drop spreading and impact are dynamic processes that also involve the interactions with solid surfaces.Therefore, fluid(air/gas or liquid)interacting with solid surfaces must also be accounted for in these simulations.This then introduces the concept of“contact line”.

    Direct simulations of fluid interacting with solid boundaries encounters a mathematical dilemma in capturing the moving contact line along a no-slip boundary:the fluid velocity is finite at the free-surface but zero on the solid surface[Dupon-t and Legendre(2010)].Therefore a stress singularity is formed.The divergent stress comes from the fact that the continuum assumption of the fluid is no longer valid within the molecular length scale near the three phase contact region.The molecular interactions between the fluid and the solid within the molecular length scale are required in the model.Besides direct simulations,there are various numerical models to deal with the contact line problem.To apply a non-slip boundary condition on a solid surface,one of the models suggests[Di and Wang(2009)]to assume there is a thin liquid film as a precursor model.The singular stress is naturally removed if the inclusion of a microscopic precursor film is placed in front of the apparent contact line.However,the assumption of the existence of the thin liquid film is under debate because if it does exist,it is generally not fast enough to stay ahead of the contact line[Chen,Mertz,and Kulenovic(2009)].Other studies have been made to address this problem by relaxing the no-slip boundary condition with a slip boundary model to resolve the singular stress when solving the Navier-Stokes equations.Due to the lack of description of the microscopic interactions between the fluid and solid in the subgrid scale,the contact angle typically needs to be specified from some empirical models to treat the contact line.Many multiphase simulating techniques such as front tracking[Muradoglu and Tasoglu(2010)],level-set[Spelt(2005);Chen,Mertz,and Kulenovic(2009)],VOF[Sikalo,Wihelm,and Roisman(2005);Afkhami,Zaleski,and Bussmann(2009)]have adopted this approach to simulate contact line problems successfully.The contact angle is either imposed to the explicit interfacial markers at the contact line(front tracking method)or specified as a boundary condition(level-set and VOF).Despite the limitation of the results been dependent on the slip-length and thus on the grid spacing,the method is still widely adopted due to its simplicity and easiness to implement into an existing multiphase flow algorithm.

    The biggest challenge here,however,is modeling thedrop impact,where a drop impacts on a solid surface and upon impact,the drop spatters into small drops.In other words,the topology,or the geometry of the interface of the drop changes as the drop breaks up.To handle topology changes numerically is not a trivial task because of the connectivities among the points that represent the interface cannot be maintained through time.The re-construction and re-connecting of the drop interface as it goes through dramatic changes can be quite daunting.In our recent development,we proposed a connectivity-free front tracking(CFFT)method[Wang,Wang,and Zhang(2013)].The main goal of this development was to handle topology changes of the interfacewithoutthe constraint of the connectivities of the interfacial points,therefore the name “connectivity-free”.It combines the merits of the front tracking and front capturing approaches such that explicit interfacial points were used for interface representation to achieve better volume conservation,and constructing a smooth and accurate indicator field based on the actual location and geometry of the interface for fluid properties(curvature and surface tension)calculations.

    Naturally,to successfully model drop spreading and impact involves two main features to be included in an algorithm:(1)a multiphase model without connectivity to handle drop breaking-up or spattering,and(2)a fluid-solid contact line model that accurately captures the contact line dynamics.The objective of this paper is to couple a dynamic contact line model into the CFFT algorithm so that drop spreading and impact can be studied and analyzed.To our knowledge,this is a first attempt in coupling a contact line model into a connectivity-free multiphase method.To demonstrate the capabilities of this coupled algorithm,several numerical experiments are examined in 2-D.Three-dimensional problems are not particularly difficult in either contact line model or CFFT algorithm individually,in fact they have been examined independently for 3-D problems[Shi,Bao,and Wang(2013);Wang,Wang,and Zhang(2013)]even though contact line model was mostly modeled with axis-symmetric and 2-D[Dupont and Legendre(2010)].However,the combination of the two involves quite complex 3-D surface reconstruction algorithm,which requires significantly more effort.Here,our main focus is to first develop the numerical framework,the 3-D implementation is to be completed in the near future.

    Under the framework of the CFFT method,the points regeneration scheme can cope with the interface topology change automatically,and construct the initial contact line by connecting points to the solid surface when the interface is close enough.The dynamic contact line model then reconstructs the contact by adjusting the explicitly represented interface.Coupling the contact line model to CFFT has many advantages:1)No connectivities required for the interfacial points which is hard to maintain when interface undergoes topology change;2)The break-up and coalescence,and the initial contact line formation can be automatically treated with a points-regeneration scheme;3)Easy implementation of the dynamic contact line model with the interface represented by explicit points;4)Avoid solving the advection equation in the VOF and level-set method which would require special treatment to deal with issues related to stability and volume conservation.

    The outline of this paper is as follows.In Section 2,the mathematical model is presented.We first review the governing equation of the connectivity-free front tracking method and then discuss the treatment of the moving contact line with proper empirical correlations between the contact angle and the capillary number.In Section 3,three numerical examples are studied to validate the algorithm.Finally,the conclusions are drawn in Section 4.

    2 Numerical Algorithm

    2.1 The connectivity-free front tracking method

    A ‘one fluid’approach is used to describe the multi- fluid flow regime[Hua,Stene,and Lin(2008)],which uses only a single set of Navier-Stokes equations with fluid properties varying across the interface without the need to treat the pressure jump condition.Both fluids(air and liquid)are considered as isothermal,incompressible flow.The surface tension force is added to the momentum equation as a singular force term.The overall governing equations are expressed as follows:

    where u is the velocity,pis the pressure,ρ and μ are the fluid density and viscosity,respectively,g is the gravity,and Fσis the surface tension force.

    The fluid properties(density ρ and viscosity μ)are calculated using an indicator defined asI=1 for liquid andI=0 for air.A linear interpolation is used to perform the calculation:

    The subscriptslandgdenote the liquid and air,respectively.

    In the front tracking method[Unverdi and Tryggvason(1992)],the indicator field near the interface is obtained by solving the Poisson’s equation:

    where Γ is the interface,Φ is an appropriate interpolation function and n is the unit normal,which requires to be obtainedaprioriutilizing the connectivity of the interfacial points.The connectivity of the interfacial points is sometimes difficult to maintain when the interface undergoes frequent topology changes,especially when the air-liquid starts to contact with a solid surface generating a three-phase contact,which requires the interface connectivity to be reconstructed constantly.It is even more complicated for 3-D cases.

    The CFFT method adopts an alternative way to cope with the necessity of connectivity[Wang,Wang,and Zhang(2013)].It still uses the explicit interfacial points to describe and update the interface where a good volume conservation can be achieved in a divergence-free incompressible flow field[Tryggvason,Bunner,Esmaeeli,Juric,Al-Rawahi,Tauber,Han,Nas,and Jan(2001)].However,by combining the basic concept of the front capturing methods such as level-set method and the VOF method,the indicator field is obtained first using the un-connected interfacial points through a special treatment,the unit normal and curvature are then calculated from the indicator field.

    The procedure of obtaining the indicator without interfacial points connectivity is as follows:obtain an approximate indicator field first and then correct this indicator field based on the current position of the interfacial points.In order to achieve this,we first define a set of ghost points()that are placed at the center of each fluid element(see Fig.1).Noting that this is different from the original CFFT where a ghost fluid mesh instead of points is constructed.Without the mesh,it adds flexibility when constructing the indicator.One can simply use more points within each fluid element to achieve better resolution.If the initial position of the interface is known,the approximate indicator for the ghost points at the beginning of the simulation is also known.

    In the front capturing method,the indicator or signed distance function is advected by solving an advection equation(6)implicitly at the current time step:

    Solely using advection equation to update indicator function may cause issues in stability and volume conservation.Here we use it to acquire the approximate indicator fieldIaat the current time stepnexplicitly by evaluating the indicator for each ghost point based on the previous time step indicator fieldIn?1and interface velocityun?1:

    Since the ghost points are defined arbitrarily and their indicator values are explicitly solved,it provides the possibility to increase or decrease the ghost point resolution when needed.

    Once the approximate indicator field,Ia,is obtained for the ghost points,we can interpolate the indicator field from the ghost points,xgto any point x,such as the fluid nodes and the interfacial points,through a proper delta function,Φ:

    Once the indicator values on the interfacial points are found,we then correct it so that the interface would have a constant indicator(level).In order to achieve this,a correction term δIp(xp)for each interfacial point,p,is added to Eq.(8)when performing the interpolation:

    Figure 1:Schematics of computational domain set up.

    where the subscriptpdenotes the interfacial points,Npis the number of interfacial points.The correction term needs to be solved for every interfacial point,p.If the indicator of all the interfacial points is set to be 0.5,i.e.I(xp)=0.5,then

    Here the subscriptp0is used to differentiate frompin the summation.Solving Eq.(10)yields the correction value δIpfor each interfacial point.Once the correction of the indicator function is solved,the final indicator field for the fixed fluid mesh can be obtained using Eq.(9).

    The surface tension force is calculated using a continuum surface tension force(CSF)approach[Brackbill,Kothe,and Zemach(1992)],which can convert the sin-gular point/surface force into volume force that spans the vicinity of the interface:

    where[I]is the jump of the indicator function and<ρ>is the average density at the interface.

    To make sure the interpolation is performed accurately for near solid surface region where the influence domain is incomplete,the reproducing kernel particle method(RKPM)interpolation scheme is chosen[Liu,Jun,and Zhang(1995);Liu,Jun,Li,Adee,and Belytschko(1995)].The use of a higher order interpolation function and the satisfaction of the reproducing condition are essential in accurately calculating the normal and curvature of the interface,hence in forming a correct contact line.The details of the RKPM implementation can be found in[Wang,Wang,and Zhang(2013)].

    Once the indicator field is constructed and the surface tension force is calculated,the Navier-Stokes equations are then solved to obtain the velocity and pressure fields.The interface is updated using the velocity at each interfacial point u(xp),which is interpolated from the fluid mesh:

    wherenis the current time step,n+1 is the next time step,and?tis the time step size.

    During the evolution of the interface,the interfacial points often need to be regenerated for the following purposes:1)to maintain sufficient points within each interfacial segment;2)to handle topology changes when both the breaking-up and coalescence occur by adding or deleting points near the contact;3)to form the initial contact line by connecting the nearby interfacial points to the solid surface.The regeneration of the interface is achieved by carefully selecting some points inside the fluid elements near the interface as the candidate points and projecting them onto the interface.The point projection scheme is described as follows.We first consider a candidate point at the position x with indicator valueI(x),then project the point onto the interface along the normal direction of the indicator.We allow x to move δx so that the indicator of the candidate point achieves 0.5.

    Figure 2:Contact line dynamics.

    Based on Eqs.(13)and(14),δx is solved and the candidate point is projected to a new position x0=x+δx.This scheme is performed several times until|I(x0)?0.5|< ε,where ε is a set tolerance,then a new interfacial point is regenerated.The interface topology changes and the formation of initial contact line are treated by selecting candidate points from alternate side of the interface to ensure the breaking-up and coalescence of the interface are properly treated.The detailed description on the points-regeneration scheme can be referred to[Wang,Wang,and Zhang(2013)].

    2.2 Dynamic contact line model

    The dynamics of the contact line in our study can be described as follows.At the microscopic length scale,the contact angle of a moving interface is a constant,which is equal to the equilibrium contact angle θe(see Fig.2).The contact angle is known for a given air-liquid-solid system.Since our numerical simulation does not directly account for the molecular interactions of the three phases at the microscopic length scale,we use the so-called ‘a(chǎn)pparent’dynamic contact angle which defines the macroscopic relation between the contact line and the solid surface.A dynamic contact line model is needed here to relate the dynamic contact angle and the microscopic interactions of the three phases.

    To apply no-slip boundary condition on the solid surface would yield a stress singularity at the contact line.To avoid this problem,a Navier-slip boundary condition is imposed so that the fluid is allowed to slide along the solid surface:

    whereuwis the fluid velocity on the solid surface,and λ is a pre-defined sliplength for a given problem.In reality,the true slip-length is in the order of the intermolecular distance,here we choose λ=0.001 for all the numerical examples,which does not allow too much movement of the contact line.

    The shape of the contact line is dependent on the contact angle formed with the surface and the intrinsic properties of the liquid.An empirical model is demanded here to reconstruct the contact line by correlating the dynamic contact angle and the capillary number.The model we use is first presented by Kistler[Kistler(1993)]and is widely used in other studies[Muradoglu and Tasoglu(2010);Mukherjee and Abraham(2007);Sikalo,Wihelm,and Roisman(2005);Roisman,Opfer,Tropea,Raessi,Mostaghimi,and Chandra(2008)].It is proven to be suitable in dealing with contact line movement:

    Cacl= μUcl/σ(Uclis the capillary number of the contact line;Uclis the velocity of the contact line. θeis the equilibrium contact angle.fHoffis the Hoffman’s function which is defined as:

    Due to hysteresis,the contact line can only advance when the contact angle is beyond the advancing contact angle θA,or recede when the contact angle is below the receding contact angle θR.If θAand θRare prescribed for a particular system,the dynamic contact angle θDcan still be computed using Eq.(16)by substituting θewith θA(when advancing)or θR(when receding).

    Under the framework of the connectivity-free front tracking method,the interface is explicitly presented by interfacial points.Once the interface gets close enough to the solid surface,the points-regeneration scheme starts to connect the interfacial points to the surface to form an initial contact line.Therefore,the remaining task is to use the empirical model to reconstruct the initial contact line and predict the correct movement.

    Figure 3:Schematics of contact line construction.

    The schematics of the contact line construction is shown in Fig.3.At each time step,if the interfacial points are found to be connected to the solid,the contact line needs to be reconstructed so that the unit normals of the interfacial points connecting to the solid surface obey the predicted dynamic contact angle using Eq.(16).To determine the dynamic contact angle θD,the contact line velocityUclmust be specified.To use the velocity of the three-phase contact point xOas the contact line velocity,an iterative procedure would be required to determine the location and the velocity of the three-phase contact point since each of them is affected by the other.The iteration would not reach convergence when the contact line velocity approaches zero[Muradoglu and Tasoglu(2010)].Here we let the tangential(to the surface)component of a reference point velocity vRto be the velocity of the contact lineUcl.The reference point xRis located at a distancedfrom the solid surface.The distance can be specified asd=ε?x,where?xis the grid spacing,ε is a scaling factor specified as 0.9 arbitrarily.This reference point is explicitly generated during the regeneration scheme.The unit normal nRand velocity vRfor this particular point can also be calculated.To determine whether the contact line is advancing or receding,simply evaluate the sign of vR·nR.The unit normal for each interfacial point is defined to be along the outward normal direction.Therefore,if vR·nR>0,then the contact line is advancing;if vR·nR<0,then the contact line is receding.

    Once the contact line velocity is determined,the dynamic contact angle,hence the unit normal for the point connected to the solid surface,can be calculated from Eq.(16).The contact line is reconstructed using an arc from a perfect circle.The advantage to use an arc to fit is that the outward unit normal can transit smoothly from the reference point xRto the end point of the contact line xO.For a 2-D construction,assuming the coordinates(xR,yR)and the unit normal(nx,ny)for the reference point are known.The unit normal of the point connected to the solid(mx,my)is calculated from the dynamic contact angle.The outward normal for the solid surface is denoted as(px,py).The goal is to find the coordinates of the last point of the contact line(xO,yO).Because the reference point and the last point are along the same circle,assuming the center of the circle is located at(xC,yC),then the following equations must be satisfied:

    To rule out any unrealistic solutions,the following condition must be satisfied:

    Based on the conditions(18)to(22),the coordinates at the point of contact,or the last point of the contact line are calculated.More interfacial points can be added between the reference point and three-phase contact point following the arc,providing a smooth reconstructed contact line.

    The method can be easily extended to the 3-D case.Whenever a reference point is selected and its unit normal is obtained,the contact line is constructed in the plane formed by the unit normal of the reference point(n)and the normal of the solid surface(p),similar to the 2-D contact construction.The actual locations of the newly constructed points can be obtained by transforming the coordinates from the 2-D plane to the 3-D space.

    Once the contact line construction is completed,a smooth indicator field is then obtained following the CFFT method.With accurate evaluations of the interface(contact line)normal by adopting the RKPM interpolation,the following condition can be automatically satisfied:

    The subscriptwpdenotes the interfacial point that directly connects the surface.nwpis the unit normal calculated from the dynamic contact angle.Condition(23)is used as a boundary condition for the front capturing method when solving the advection equation.Our approach achieves this condition using a different manner,as described in this section.

    2.3 Numerical procedure

    The overall numerical procedure can be summarized as follows.

    1)Perform the points-regeneration scheme to regulate the interfacial points,deal with the topology changes such as breaking-up and coalescence if occur,and connect the interface to the solid surface when interface gets sufficiently close.

    2)If a contact line is detected,calculate the contact angle θ.If θR< θ < θA,the contact line movement is under hysteresis,there is no need to reconstruct.If θ > θAor θ < θR,reconstruct the contact line following the procedure mentioned above.3)Construct a smooth indicator field based on the interface using the CFFT method.4)Calculate the fluid properties and the surface tension force using Eqs.(3),(4)and(11).Solve the Navier-Stokes equations to obtain the velocity and pressure fields for the fluid.

    5)Finally,update the interface based on the velocity field asThen go back to step 1)for the next time step.

    3 Numerical Examples

    The examples to be studied here use non-dimensionalized dimensions and material parameters.

    3.1 Drop spreads on a horizontal surface

    This numerical study is to examine the spreading of a drop on a horizontal solid surface due to gravity.Opposing to the spreading,the surface tension force tries to maintain the sphericity of the drop.The movement of the contact line is constrained by the dynamic contact angle,where the drop eventually reaches an equilibrium shape.With a given static or equilibrium contact angle(θe),the final shape of the drop can be uniquely determined.

    A 2-D circular drop is initially placed on a horizontal surface.The schematic is shown in Fig.4.The computational domain is a 6×2 rectangular box,which is discretized using 53,200 uniform quadrilateral elements.The drop with a radius of 0.5 is placed at the bottom center of the box.The density and viscosity ratios of the liquid drop and the ambient air are ρl/ρg=1000 and μl/μg=100,respectively.To analyze the final shape of the drop,we perform the simulations with and without the dynamic contact line model.For the cases with dynamic contact line model,we perform the simulation for 3 different equilibrium contact angles(θe=60?,90?and 120?)with Eotvos numberEo=1.0.The Eotvos number,Eo= ρlg/σ,represents the ratio of gravity and capillary force,where ρlis the density of the drop,R0is the radius of the initial drop,gis the gravity,and σ is the surface tension.

    Figure 4:Schematic of a drop spreading on a horizontal plane.

    In this example,surface hysteresis is not considered.The final shapes of the drop at equilibrium of these 4 cases (one without the dynamic contact line model and 3 with the dynamic contact line model)are shown in Fig.5.We can observe that when the dynamic contact line model is used,the final shapes of the drop are successfully reproduced for their corresponding equilibrium contact angles.However,when the dynamic contact line model is not used,the drop shape is purely controlled by the gravity and the surface tension force.Without any mechanism to constrain the contact line,the drop shape does not account for any wettability of the surface.In this case,the equilibrium contact angle converges to around 100?.

    In order to analyze the effect ofEonumber on the final shapes of the spreading drop,the simulations are performed forEonumber varying from 0.01 to 50,which is achieved by adjusting the ratio of gravity and surface tension.

    To further analyze the effect of theEonumber on the final shapes of the spreading drop,the simulations are performed forEonumber varying from 0.01 to 50,which is achieved by adjusting the ratio of gravity and surface tension.WhenEo<<1,the capillary force has a dominant role in the movement of the drop,the shape of the drop cap remains round.According to[Dupont and Legendre(2010)],the maximum height of the dropH0measured from the solid surface is independent ofEo,and can be determined as a function of the equilibrium contact angle as:

    WhenEo>>1,the gravity controls the shape of the drop and the maximum heightH∞of the drop is a function of both θeandEonumber:

    Figure 5:Final shapes of the drop for different equilibrium contact angles and without dynamic contact line model.

    Figure 6:Normalized maximum height vs.Eo.

    Figure 7:Schematic of a drop moving on an oblique surface.

    Here we compare the normalized maximum height of the drop(H/H0)at equilibrium for a large range of Eotvos number(0.01,0.05,0.1,0.25,0.5,1.0,2.5,5.0,10,50)on a hydrophilic(θe=60?)and a hydrophobic(θe=120?)solid surfaces with analytical prediction(Eqs.(24)and(25))in Fig.(6).It can be clearly seen that the computed drop height agrees very well with the asymptotic solutions obtained from Eqs.(24)and(25).The transition of the final shape of the drop from a round cap to a puddled shape asEonumber increased is clearly seen.This simple test shows that the dynamic contact line model that is coupled into the CFFT algorithm can accurately construct the contact line according to a given contact angle.The final shape of the spreading drop can be successfully obtained.

    3.2 Drop spreads on an oblique surface

    In this example,a drop with hysteresis moving on an oblique surface with variousEonumber is simulated.When the surface inclination angle is at a critical inclination angle αC,the gravity and capillary force balance with each other,which is also called the critical condition(See Fig.7).If the inclination angle of the surface is larger than αC,both the front and the trailing edges of the drop would start to move.The critical inclination angle is also influenced by the given hysteresis angles,and can be calculated as[Dupont and Legendre(2010)]:

    whereAis the entire surface area of the 2-D drop. θAand θRare the receding and advancing angles,respectively.The upper caseAandRdenote the predefined advancing and receding contact angles based on hysteresis.We use the lower caseaandrto describe the actual advancing and receding contact angles obtained from the numerical calculation based on the dynamic contact line model.Re-arranging the equation withEonumber we get:

    In this set-up,a semi-circular drop is placed on a solid surface.The schematic is the same as previous numerical example.The inclination angle of the solid surface is adjusted by manipulating the direction of the gravity.The density and viscosity ratios of the liquid drop and air are 1000 and 100,respectively.Three sets ofEonumbers 0.5,1 and 2 and two sets of hysteresis angles(θA,θR)=(100?,80?)and(120?,60?)are chosen.The equilibrium contact angle is set to be 90?.For each combination of anEonumber and a set of hysteresis angles,a series of inclination angles with an increment of 5?or 10?near and away from the estimated critical inclination angle are examined using the simulation model.Each simulation is run until the system reaches equilibrium.From the simulations,a critical inclination angle is identified by linear interpolation of inclination angles right before and after the the drop reaches hysteresis.An example is shown in Fig.8 for the combination

    Eo=1 and(θA,θR)=(120?,60?).This figure shows the variation of the advancing and receding contact angles for a range of selected inclination angles α from 10?to 90?.When the inclination angle increases,the drop starts to deform in order to maintain a balance between the gravity and the capillary force.The advancing contact angle reaches beyond the hysteresis angle(θA=120?) first and exceeds a little during the rest of the simulation,which agrees with the experimental observation in[ElSherbini and Jacobi(2006)].Similar numerical simulations[Dupont and Legendre(2010)]also confirm this conclusion.The receding contact line starts to move around α =39?,at which point the drop starts to slide due to both the advancing and receding contact lines are moving.The drop shapes when α =20?,40?and 80?are shown in Fig.9.The evolution of the drop and the sliding behavior after the drop reaches the critical condition(α =80?)can be clearly seen.

    Figure 8:Advancing and receding contact angle for different inclination angles for Eo=1 and(θA,θR)=(120?,60?).

    Figure 9:Drop shapes at α =20?,40? and 80? for Eo=1 and(θA,θR)=(120?,60?).

    Figure 10:sin αCvs.2(cosθR?cos θA)/πEo.

    Table 1:Calculated critical angle.

    Figure 11:Normalized receding contact angle vs.Bo before the drop reaches critical condition.The experimental data fit is obtained from[ElSherbini and Jacobi(2006)].

    The critical angle identified from the simulation for differentEonumber and hysteresis angles are listed in Table(1).We also compare the numerical and theoretical critical angles obtained from Eq.(27),shown in Fig.10.A very good agreement is achieved.For the combination ofEo=0.5 and(θA,θR)=(120?,60?),there is no motion for the drop.That is because during the entire simulation the receding contact angle θrnever reaches the predefined receding contact angle of θR=60?,and the capillary force completely balances with the gravity.Therefore the drop hangs onto the surface.

    In a previous work,ElSherbini and Jacobi[ElSherbini and Jacobi(2006)]found a linear relationship between the minimum receding contact angle of a 3-D liquid drop andBonumberbefore the drop reaches the critical condition:

    Figure 12:Schematics of a drop impacts on a solid surface.

    This relationship is based on the experimental data using various testing surface,liquid,drop size,and inclination angle combinations,which result in different hysteresis angles and liquid properties.To see if our simulation can reproduce this relationship,we plot the non-dimensionalized receding contact angle θr/θAfor variousBonumbers before the drop reaches the critical condition in Fig.11.We may notice that θr/θAdecreases almost linearly whenBoincreases,which is consistent with the linear relationship obtained from the experimental data.However,the experiment is performed for 3-D drops with θrbeen the minimum receding contact angle and the experiment θAbeen restricted to from 49?to 112?.Those factors could be the reasons for the small discrepancy of the simulated results and the linear fit from the experimental data.

    3.3 Drop impacts on a solid

    To show the robustness of the implementation,we simulate a more complicated case with a drop impacting onto a round solid obstacle.The geometry setup is shown in Fig.12.The computational domain is a 6×12 rectangular box,which is discretized with 11,1861 quadrilateral elements.A round obstacle is placed at the center of the box with a diameter of 2.0.A drop with a diameter ofD=1.6 is placed at a distance of 10.0 measured from the bottom.All the boundaries are set as Navier-slip boundary conditions.The density and viscosity ratios of the drop and the ambient air are 1000 and 100,respectively.The time step size is 0.001.The gravity acts along the negative vertical direction.The shape of the drop when moving in the ambient air is governed by the gravity and surface tension force.The non-dimensional Reynolds andEonumbers are defined as:

    Figure 13:Shapes of the drop when impacting on the solid with Re=14.1,Eo=12.8.

    Figure 14:Shapes of the drop when impaction on the solid with Re=141,Eo=12.8.

    The values of gravity and surface tension force are manipulated so that two Reynolds number cases withEo=12.8 are studied:1)Re=14.1;2)Re=141.The equilibrium contact angle for both cases is set to be 90?.The advancing and receding angles(θA,θR)are set as(100?,80?),respectively.

    Figs.13 and 14 show the drop shapes at different non-dimensional times for case 1 and case 2,respectively.Due to gravity,the drop falls and impacts onto a round obstacle.When the drop is close enough to the solid surface,the points-regeneration scheme helps form the contact line by connecting the interfacial points to the solid obstacle(Figs.13(b)and 14(b)).Then the movement of the contact line conforms to the dynamic contact line model(using Eq.(16))and the given surface hysteresis.For case 1 that has a relatively small Reynolds number,the drop tends to drip down(Fig.13(c)).When the surface tension force can no longer be balanced with the gravity of the drips,the drop breaks up into small ones.As shown in Fig.13,the drop undergoes two stages of breaking-up and forms two large drops and two small ones(Fig.13(d)).The two large drops coalesce during falling(Fig.13(f))and eventually impact on the bottom surface and form a round cap(Fig.13(h)).Due to the accumulative precision error,the small drops are not formed at the exact same time.Both of them eventually merge with the two large drops.In case 2 where the Reynolds number is 141,as shown in Fig.14,a different scenario happens.The original drop hits the round obstacle with a larger impact velocity.Therefore the impact is much more dramatic compared to the low Reynolds number case.The drop breaks up into small ones that slash around.Instead of forming a round cap at the bottom surface,the topology of the interface changes constantly.It even climbs onto the side walls during the spreading of these broken drops.

    4 Conclusions

    In this work,we propose and set up a numerical framework to model drop spreading and impact by coupling a connectivity-free front tracking(CFFT)method with a dynamic contact line model with surface hysteresis.This numerical method can be used as a tool to perform bloodstain pattern analysis as it aims at accurately predicting the drop behaviors such as motion,impact and spreading under gravity,viscosity and surface tension force.The CFFT method explicitly represents the interface using interfacial points without the need to maintain the connectivities.Therefore,the topology change can be easily handled for blood drop spattering.In order to provide more accurate simulation for contact line problem,we adopt a macroscopic dynamic contact line model to predict the contact angle.This model is implemented by reconstructing the contact region of the air-liquid interface to adjust the actual geometrical contact angle according to the model prediction.With the explicit interfacial points representation,the reconstruction is accurately performed and the dynamic contact angle constraint is successfully imposed.

    Several test cases are shown to validate the method.In the drop spreading due to gravity case,the normalized maximum height of the drop agrees very well with the theoretical prediction.By performing the drop movement on an oblique surface case,we find the critical angles of the oblique surface with variousEonumber and surface hysteresis to match well with the theoretical ones,and numerical results of the receding angles for the inclined drop agree with the experimental data.Finally,the method is applied to a more complicated case by impacting the drop onto a solid obstacle.The drop undergoes frequent interface topology changes.The formation and movement of the contact lines when the drop interacts with solid surfaces are well presented.The results show the different behaviors of the drop due to different Reynolds numbers hence the impact velocities,which indicates great flexibility and robustness of the method.

    We acknowledge that the implementation as well as the test cases were limited to 2-D at this point due to the complex and heavy computational cost involved in constructing 3-D surfaces without interfacial connectivity.However,this work is to demonstrate the applicability in accurately modeling drop spreading and impact problems that are not trivial in many numerical aspects.To our knowledge,this is a first attempt in combining a connectivity-free multiphase approach with a dynamic contact line model,which shows great promises in building a robust numerical framework to be widely adopted in BPA and forensic science.

    Afkhami,S.;Zaleski,S.;Bussmann,M.(2009):A mesh-dependent model for applying dynamic contact angles to vof simulations.Journal of Computational Physics,vol.207,pp.5370–5389.

    Attinger,D.;Moore,C.;Donaldson,A.;Jafari,A.;Stone,H.A.(2013):Fluid dynamics topics in bloodstain pattern analysis:Comparative review and research opportunities.Forensic Science International,vol.231,no.1-3,pp.375–396.

    Brackbill,J.U.;Kothe,D.B.;Zemach,C.(1992): A continuum method for modeling surface tension.Journal of Computational Physics,vol.100,pp.335–354.

    Cervone,A.;Manservisi,S.;Scardovelli,R.(2010):Simulation of axisymmetric jets with a finite element Navier-Stokes solver and a multilevel VOF approach.Journal of Computational Physics,vol.229,no.19,pp.6853–6873.

    Chen,Y.;Mertz,R.;Kulenovic,R.(2009): Numerical simulation of bubble formation on orifice plates with moving contact line.International Journal of Multiphase Flow,vol.35,pp.66–77.

    Di,Y.;Wang,X.(2009):Precursor simulations in spreading using a multi-mesh adaptive finite element method.Journal of Computational Physics,vol.228,pp.1380–1390.

    Dupont,J.-B.;Legendre,D.(2010):Numerical simulation of static and slideing drop with contact angle hysteresis.Journal of Computational Physics,vol.229,pp.2453–2478.

    ElSherbini,A.;Jacobi,A.(2006):Retention forces and contact angles for critical liquid drops on non-horizontal surfaces.Journal of Colloid and Interface Science,vol.299,pp.841–849.

    Herrmann,M.(2008): A balanced force refined level set grid method for twophase flows on unstructured flow solver grids.Journal of Computational Physics,vol.227,no.4,pp.2674–2706.

    Hirt,C.;Nichols,B.(1981):Volume of fluid(vof)method for the dynamics of free boundaries.Journal of Computational Physics,vol.39,no.1,pp.201–225.

    Hua,J.;Stene,J.F.;Lin,P.(2008):Numerical simulation of 3d bubbles rising in viscous liquids using a front tracking method.Journal of Computational Physics,vol.227,pp.3358–3382.

    Karger,B.;Rand,S.;Fracasso,T.;Pfeiffer,H.(2008): Bloodstain pattern analysis-casework experience.Forensic Science International,vol.181,no.1-3,pp.15–20.

    Kistler,S.(1993): Hydrodynamics of wetting.In Berg,J.(Ed):Wettability,pp.311–429,New York.Marcel Dekker.

    Liu,W.K.;Jun,S.;Li,S.;Adee,J.;Belytschko,T.(1995):Reproducing kernel particle methods for structural dynamics.International Journal for Numerical Methods in Engineering,vol.38,pp.1655–1679.

    Liu,W.K.;Jun,S.;Zhang,Y.F.(1995):Reproducing kernel particle methods.International Journal for Numerical Methods in Fluids,vol.20,pp.1081–1106.

    Mukherjee,S.;Abraham,J.(2007):Investigations of drop impact on dry walls with a lattice-baltzmann model.Journal of Colloid and Interface Science,vol.312,pp.341–354.

    Muradoglu,M.;Tasoglu,S.(2010):A front-tracking method for computational modeling of impact and spreading of viscous droplets on solid walls.Computers&Fluids,vol.39,pp.615–625.

    Muradoglu,M.;Tryggvason,G.(2008): A front-tracking method for computation of interfacial flows with soluble surfactants.Journal of Computational Physics,vol.227,no.4,pp.2238–2262.

    Roisman,I.;Opfer,L.;Tropea,C.;Raessi,M.;Mostaghimi,J.;Chandra,S.(2008): Drop impact onto a dry surface:Role of the dynamic contact angle.Colloids and Surfaces A:Physicochemical and Engineering Aspects,vol.322,pp.183–191.

    Shi,Y.;Bao,K.;Wang,X.-P.(2013): 3d adaptive finite element method for a phase field model for the moving contact line problems.Inverse Problems and Imaging,vol.7,no.3,pp.947–959.

    Sikalo,S.;Wihelm,H.-D.;Roisman,I.(2005): Dynamic contact angle of spreading droplets:Experiments and simulations.Physics of Fluids,vol.17,pp.615–625.

    Spelt,P.D.(2005): A level-set approach for simulations of flows with multiple moving contact lines with hysteresis.Journal of Computational Physics,vol.207,pp.389–404.

    Tryggvason,G.;Bunner,B.;Esmaeeli,A.;Juric,D.;Al-Rawahi,N.;Tauber,W.;Han,J.;Nas,S.;Jan,Y.-J.(2001):A front-tracking method for the computations of multiphase flow.Journal of Computational Physics,vol.169,no.2,pp.708–759.

    Unverdi,S.O.;Tryggvason,G.(1992): A front-tracking method for viscous,incompressible,multi- fluid flows.Journal of Computational Physics,vol.10,pp.25–37.

    Wang,C.;Wang,X.;Zhang,L.(2013):Connectivity-free front tracking method for multiphase flows with free surfaces.Journal of Computational Physics,vol.241,pp.58–75.

    Witteveen,J.;Koren,B.;Bakker,P.(2007):An improved front tracking method for the euler equations.Journal of Computational Physics,vol.224,no.2,pp.712–728.

    Xu,J.-J.;Li,Z.;Lowengrub,J.;Zhao,H.(2006): A level-set method for interfacial flows with surfactant.Journal of Computational Physics,vol.212,no.2,pp.590–616.

    亚洲国产日韩欧美精品在线观看 | 色av中文字幕| 欧美日韩亚洲国产一区二区在线观看| 久久国产精品人妻蜜桃| 国产精品99久久99久久久不卡| ponron亚洲| 日韩欧美三级三区| 啪啪无遮挡十八禁网站| 中出人妻视频一区二区| 99国产精品99久久久久| 12—13女人毛片做爰片一| 在线观看一区二区三区| 亚洲精品中文字幕一二三四区| 亚洲成av人片免费观看| 高清在线国产一区| 国产黄片美女视频| 美女午夜性视频免费| 一个人免费在线观看的高清视频| 午夜福利视频1000在线观看| 精品久久久久久久久久久久久| netflix在线观看网站| 麻豆成人av在线观看| 国产野战对白在线观看| 国产精品自产拍在线观看55亚洲| 欧美成狂野欧美在线观看| svipshipincom国产片| 亚洲自偷自拍图片 自拍| 日韩精品青青久久久久久| 亚洲国产欧洲综合997久久,| 毛片女人毛片| 窝窝影院91人妻| 日本黄色视频三级网站网址| 欧美午夜高清在线| 免费搜索国产男女视频| 精品无人区乱码1区二区| 怎么达到女性高潮| 少妇粗大呻吟视频| 九色成人免费人妻av| 午夜精品在线福利| 国产在线精品亚洲第一网站| 中文亚洲av片在线观看爽| 国内揄拍国产精品人妻在线| 亚洲黑人精品在线| 国产激情久久老熟女| 少妇熟女aⅴ在线视频| 亚洲成人精品中文字幕电影| 欧美日韩亚洲综合一区二区三区_| 给我免费播放毛片高清在线观看| 午夜福利成人在线免费观看| 又黄又粗又硬又大视频| 亚洲自拍偷在线| av在线天堂中文字幕| 美女 人体艺术 gogo| 中文字幕久久专区| 悠悠久久av| 美女黄网站色视频| 18禁裸乳无遮挡免费网站照片| 国产单亲对白刺激| 又黄又粗又硬又大视频| 18禁黄网站禁片免费观看直播| 亚洲精品美女久久av网站| 香蕉国产在线看| 天堂影院成人在线观看| 欧美日韩精品网址| 久久这里只有精品中国| 国产精品久久久久久亚洲av鲁大| 亚洲精品美女久久av网站| 亚洲精品久久国产高清桃花| 久久午夜综合久久蜜桃| 欧美日韩黄片免| 国产高清视频在线观看网站| 给我免费播放毛片高清在线观看| 级片在线观看| 精品国产超薄肉色丝袜足j| 国产精品久久久久久精品电影| 成人午夜高清在线视频| 蜜桃久久精品国产亚洲av| 久久久久国产精品人妻aⅴ院| 日韩欧美国产一区二区入口| 亚洲国产欧美人成| 精品一区二区三区四区五区乱码| 免费观看精品视频网站| 国产爱豆传媒在线观看 | 精品国内亚洲2022精品成人| 伊人久久大香线蕉亚洲五| 变态另类丝袜制服| 国产成人一区二区三区免费视频网站| 两个人看的免费小视频| 色哟哟哟哟哟哟| 88av欧美| 久久香蕉激情| 成人精品一区二区免费| 人人妻,人人澡人人爽秒播| 久久久久久久午夜电影| 很黄的视频免费| 欧美最黄视频在线播放免费| 嫩草影视91久久| a级毛片a级免费在线| 亚洲人成伊人成综合网2020| 无限看片的www在线观看| 日韩成人在线观看一区二区三区| a级毛片a级免费在线| 在线看三级毛片| 中文资源天堂在线| 97人妻精品一区二区三区麻豆| 九九热线精品视视频播放| aaaaa片日本免费| 又黄又粗又硬又大视频| 国产av一区二区精品久久| 久久精品成人免费网站| 91大片在线观看| 精品国产超薄肉色丝袜足j| 九色成人免费人妻av| 午夜激情福利司机影院| 又黄又粗又硬又大视频| 99久久精品热视频| 真人一进一出gif抽搐免费| 国内精品久久久久精免费| 黄色 视频免费看| 床上黄色一级片| 99久久国产精品久久久| 欧美绝顶高潮抽搐喷水| 黄色女人牲交| 午夜精品在线福利| 久久久久久久久中文| 2021天堂中文幕一二区在线观| 亚洲狠狠婷婷综合久久图片| 成年女人毛片免费观看观看9| 香蕉久久夜色| 男人舔女人下体高潮全视频| 久久性视频一级片| 亚洲av片天天在线观看| 男人舔女人下体高潮全视频| 在线观看免费日韩欧美大片| 看片在线看免费视频| 12—13女人毛片做爰片一| 国产精品,欧美在线| 一边摸一边做爽爽视频免费| 身体一侧抽搐| av在线播放免费不卡| 黄色女人牲交| 男女之事视频高清在线观看| 日本一区二区免费在线视频| 看免费av毛片| 国产伦在线观看视频一区| 国产高清视频在线播放一区| 美女午夜性视频免费| 九九热线精品视视频播放| 桃红色精品国产亚洲av| 国产v大片淫在线免费观看| 日韩大尺度精品在线看网址| 色老头精品视频在线观看| 嫩草影视91久久| 他把我摸到了高潮在线观看| 国产精品爽爽va在线观看网站| 国产一区二区在线av高清观看| 久久天堂一区二区三区四区| 久久久久久九九精品二区国产 | 成人高潮视频无遮挡免费网站| 欧美性长视频在线观看| 久久久久久免费高清国产稀缺| 99久久久亚洲精品蜜臀av| 在线视频色国产色| 亚洲国产欧美人成| 日韩大码丰满熟妇| 在线看三级毛片| 免费在线观看日本一区| 午夜激情av网站| 久久国产乱子伦精品免费另类| 精品乱码久久久久久99久播| 黄频高清免费视频| 可以在线观看的亚洲视频| 亚洲欧美日韩高清在线视频| 国产精品一区二区精品视频观看| 久久精品人妻少妇| 色综合站精品国产| 99精品久久久久人妻精品| 久久人妻av系列| 久久婷婷成人综合色麻豆| 国产精品精品国产色婷婷| 欧美成人免费av一区二区三区| 欧美绝顶高潮抽搐喷水| 国产精品免费一区二区三区在线| 久久九九热精品免费| 99久久99久久久精品蜜桃| 99re在线观看精品视频| 国产真实乱freesex| 午夜日韩欧美国产| 他把我摸到了高潮在线观看| xxxwww97欧美| 欧美成人午夜精品| 亚洲性夜色夜夜综合| 国产成人啪精品午夜网站| 身体一侧抽搐| 久久中文字幕人妻熟女| 黄色丝袜av网址大全| 老熟妇仑乱视频hdxx| 国产黄色小视频在线观看| 九色国产91popny在线| 好男人电影高清在线观看| 黄频高清免费视频| 一级毛片女人18水好多| 五月伊人婷婷丁香| 成人手机av| 女人爽到高潮嗷嗷叫在线视频| 村上凉子中文字幕在线| 午夜精品在线福利| 亚洲国产精品合色在线| 深夜精品福利| 少妇的丰满在线观看| 亚洲人成网站在线播放欧美日韩| 国产一区二区激情短视频| 狂野欧美激情性xxxx| 日本成人三级电影网站| 最近视频中文字幕2019在线8| 91麻豆av在线| 每晚都被弄得嗷嗷叫到高潮| 欧美成人性av电影在线观看| 一级毛片女人18水好多| 亚洲成人中文字幕在线播放| 国产精华一区二区三区| 又黄又粗又硬又大视频| 国产欧美日韩精品亚洲av| 首页视频小说图片口味搜索| 亚洲国产高清在线一区二区三| 床上黄色一级片| 精品久久蜜臀av无| 午夜视频精品福利| 国产精品久久久久久人妻精品电影| 国产成人啪精品午夜网站| 悠悠久久av| 18禁美女被吸乳视频| 久久久精品大字幕| 国产一区二区在线观看日韩 | 精品国产亚洲在线| 1024香蕉在线观看| 亚洲五月婷婷丁香| 国产探花在线观看一区二区| 亚洲自拍偷在线| 男人的好看免费观看在线视频 | 国产一区二区三区在线臀色熟女| 一卡2卡三卡四卡精品乱码亚洲| 亚洲专区国产一区二区| 免费观看人在逋| 成人18禁在线播放| 2021天堂中文幕一二区在线观| 国产精品久久久久久人妻精品电影| 50天的宝宝边吃奶边哭怎么回事| 欧美又色又爽又黄视频| 亚洲一码二码三码区别大吗| 久久 成人 亚洲| 在线观看一区二区三区| 非洲黑人性xxxx精品又粗又长| 国产精品国产高清国产av| 国产精品,欧美在线| 欧美丝袜亚洲另类 | 精品久久久久久,| 99久久无色码亚洲精品果冻| 国产精品永久免费网站| 日韩欧美国产一区二区入口| 色播亚洲综合网| 亚洲精品国产一区二区精华液| 动漫黄色视频在线观看| 国产精品 国内视频| 亚洲av片天天在线观看| 欧美+亚洲+日韩+国产| 欧美日韩乱码在线| 国产激情欧美一区二区| 午夜免费观看网址| 熟女少妇亚洲综合色aaa.| 日韩欧美国产一区二区入口| 国产欧美日韩一区二区精品| av福利片在线观看| 少妇裸体淫交视频免费看高清 | 丁香欧美五月| 99热6这里只有精品| 国产精品久久久久久久电影 | 久久精品影院6| 亚洲男人的天堂狠狠| 亚洲av成人不卡在线观看播放网| av超薄肉色丝袜交足视频| 亚洲av成人一区二区三| 日韩三级视频一区二区三区| 最新美女视频免费是黄的| 精华霜和精华液先用哪个| 操出白浆在线播放| 动漫黄色视频在线观看| 18美女黄网站色大片免费观看| 国产精品电影一区二区三区| 国产91精品成人一区二区三区| 非洲黑人性xxxx精品又粗又长| 国产精品99久久99久久久不卡| 亚洲精品美女久久av网站| 国产精品 欧美亚洲| 亚洲成人精品中文字幕电影| 热99re8久久精品国产| 50天的宝宝边吃奶边哭怎么回事| 欧美+亚洲+日韩+国产| 香蕉av资源在线| 精品午夜福利视频在线观看一区| 夜夜看夜夜爽夜夜摸| a级毛片a级免费在线| 精品久久久久久久久久久久久| 亚洲 国产 在线| 午夜免费激情av| 99国产精品一区二区三区| 久久香蕉国产精品| 中文字幕久久专区| 日本精品一区二区三区蜜桃| 级片在线观看| 在线永久观看黄色视频| svipshipincom国产片| 亚洲欧美日韩无卡精品| av片东京热男人的天堂| 国产精品久久久久久亚洲av鲁大| 一本一本综合久久| 黑人巨大精品欧美一区二区mp4| 亚洲人成77777在线视频| 人妻丰满熟妇av一区二区三区| 日韩精品免费视频一区二区三区| 99久久国产精品久久久| 大型黄色视频在线免费观看| 国产成人aa在线观看| 动漫黄色视频在线观看| 国产精品精品国产色婷婷| 国产99久久九九免费精品| 精品国产美女av久久久久小说| 久久久久性生活片| 伊人久久大香线蕉亚洲五| 午夜精品在线福利| 日本黄大片高清| 久久久久九九精品影院| 亚洲中文字幕日韩| 亚洲精华国产精华精| 久久久久国产精品人妻aⅴ院| 亚洲中文av在线| 国产精品一区二区免费欧美| 18禁国产床啪视频网站| 一本精品99久久精品77| 久久天堂一区二区三区四区| 欧美大码av| 国产蜜桃级精品一区二区三区| 国产69精品久久久久777片 | 国产亚洲欧美在线一区二区| 亚洲 欧美一区二区三区| 国产91精品成人一区二区三区| www国产在线视频色| 波多野结衣高清无吗| 后天国语完整版免费观看| 亚洲人成网站在线播放欧美日韩| 日本a在线网址| 国产亚洲精品一区二区www| 色综合站精品国产| 精品高清国产在线一区| 精品久久久久久久毛片微露脸| 精品第一国产精品| 岛国视频午夜一区免费看| 最近最新中文字幕大全电影3| 欧美在线黄色| 成年女人毛片免费观看观看9| 久久精品人妻少妇| 中文字幕精品亚洲无线码一区| 亚洲国产高清在线一区二区三| 51午夜福利影视在线观看| 中文资源天堂在线| 国产精品爽爽va在线观看网站| 91大片在线观看| 亚洲激情在线av| 法律面前人人平等表现在哪些方面| 2021天堂中文幕一二区在线观| 成年人黄色毛片网站| 精品国产超薄肉色丝袜足j| 一本久久中文字幕| 国产一级毛片七仙女欲春2| 给我免费播放毛片高清在线观看| av免费在线观看网站| 免费搜索国产男女视频| 精品少妇一区二区三区视频日本电影| 久久久久久久精品吃奶| 免费无遮挡裸体视频| 亚洲av熟女| 少妇熟女aⅴ在线视频| 久久精品国产亚洲av高清一级| 国产成人精品久久二区二区免费| 中文字幕人妻丝袜一区二区| 亚洲 欧美一区二区三区| 长腿黑丝高跟| 中文字幕av在线有码专区| 亚洲无线在线观看| 欧美乱码精品一区二区三区| 制服丝袜大香蕉在线| 免费在线观看视频国产中文字幕亚洲| 亚洲性夜色夜夜综合| 亚洲成av人片在线播放无| 90打野战视频偷拍视频| 91九色精品人成在线观看| 亚洲av电影在线进入| 天堂av国产一区二区熟女人妻 | www国产在线视频色| 亚洲国产精品sss在线观看| 特大巨黑吊av在线直播| 日本 av在线| 国产精品 欧美亚洲| 久久99热这里只有精品18| 欧美国产日韩亚洲一区| 日韩欧美三级三区| 日日摸夜夜添夜夜添小说| 国产一区二区激情短视频| 成人精品一区二区免费| 国内毛片毛片毛片毛片毛片| 成人国产综合亚洲| 亚洲精品久久国产高清桃花| 老汉色∧v一级毛片| 好男人在线观看高清免费视频| 亚洲人成网站在线播放欧美日韩| www.www免费av| 中文亚洲av片在线观看爽| 两性夫妻黄色片| 久久婷婷成人综合色麻豆| 欧美日韩福利视频一区二区| 国产精品影院久久| 亚洲欧美精品综合久久99| 亚洲美女视频黄频| 精品一区二区三区av网在线观看| 国产高清视频在线观看网站| 国产成人精品久久二区二区免费| 露出奶头的视频| 国产亚洲精品一区二区www| 搡老妇女老女人老熟妇| 在线观看66精品国产| 又大又爽又粗| 两个人的视频大全免费| 日本撒尿小便嘘嘘汇集6| 国产成人系列免费观看| √禁漫天堂资源中文www| 免费电影在线观看免费观看| 黄色毛片三级朝国网站| 久久久国产精品麻豆| 国产伦在线观看视频一区| 成人午夜高清在线视频| 亚洲熟女毛片儿| 午夜精品久久久久久毛片777| 精品国产亚洲在线| 久久精品夜夜夜夜夜久久蜜豆 | 人妻夜夜爽99麻豆av| 午夜精品在线福利| 波多野结衣高清无吗| 久久精品国产99精品国产亚洲性色| 日韩欧美国产一区二区入口| 首页视频小说图片口味搜索| 丰满的人妻完整版| 黄色毛片三级朝国网站| 搡老妇女老女人老熟妇| 午夜两性在线视频| 非洲黑人性xxxx精品又粗又长| 亚洲精品在线观看二区| 亚洲国产中文字幕在线视频| 亚洲男人天堂网一区| 国产1区2区3区精品| 极品教师在线免费播放| 老司机在亚洲福利影院| 久久午夜亚洲精品久久| 婷婷精品国产亚洲av在线| 桃色一区二区三区在线观看| 少妇被粗大的猛进出69影院| 日日夜夜操网爽| 丝袜美腿诱惑在线| 在线观看66精品国产| 午夜精品久久久久久毛片777| 成人国产一区最新在线观看| 欧美极品一区二区三区四区| 免费在线观看完整版高清| 精品久久久久久久人妻蜜臀av| 极品教师在线免费播放| 亚洲专区国产一区二区| 日本一区二区免费在线视频| www日本黄色视频网| 精品久久久久久久久久久久久| 99精品欧美一区二区三区四区| 免费在线观看日本一区| 一夜夜www| 一级黄色大片毛片| 亚洲成人久久爱视频| 亚洲成人中文字幕在线播放| 99热这里只有是精品50| av在线天堂中文字幕| 国产精品1区2区在线观看.| 一本久久中文字幕| 欧美又色又爽又黄视频| 国产精品 欧美亚洲| 亚洲精品美女久久av网站| 国产成+人综合+亚洲专区| 久久精品成人免费网站| 男人舔奶头视频| 成人三级做爰电影| 老汉色av国产亚洲站长工具| 波多野结衣巨乳人妻| 午夜福利免费观看在线| 中文字幕高清在线视频| 久久久久久大精品| 午夜日韩欧美国产| 91av网站免费观看| 日韩国内少妇激情av| 在线观看美女被高潮喷水网站 | 欧美精品啪啪一区二区三区| 日韩av在线大香蕉| 老熟妇乱子伦视频在线观看| 国产精品 国内视频| 777久久人妻少妇嫩草av网站| av视频在线观看入口| 亚洲成av人片在线播放无| 一二三四社区在线视频社区8| 精品午夜福利视频在线观看一区| 最好的美女福利视频网| 最新在线观看一区二区三区| 2021天堂中文幕一二区在线观| 一进一出抽搐动态| 精品电影一区二区在线| 在线观看日韩欧美| 好男人电影高清在线观看| 国产精品美女特级片免费视频播放器 | 国产亚洲av嫩草精品影院| 在线观看免费视频日本深夜| 丝袜人妻中文字幕| 亚洲成人免费电影在线观看| 国产精品久久久久久人妻精品电影| 一级毛片高清免费大全| 91麻豆精品激情在线观看国产| 亚洲免费av在线视频| 男人舔女人的私密视频| 黄色视频,在线免费观看| 可以在线观看毛片的网站| 女生性感内裤真人,穿戴方法视频| 国内精品久久久久精免费| 亚洲熟妇熟女久久| 天堂影院成人在线观看| 少妇被粗大的猛进出69影院| 亚洲欧美精品综合久久99| 日韩欧美一区二区三区在线观看| 97超级碰碰碰精品色视频在线观看| 国产精品影院久久| 国产精华一区二区三区| 午夜激情av网站| 久久精品91无色码中文字幕| 日本熟妇午夜| 国产成人精品久久二区二区91| 一级毛片女人18水好多| 久久久久久国产a免费观看| 男女那种视频在线观看| 九色成人免费人妻av| www国产在线视频色| 亚洲精品在线观看二区| cao死你这个sao货| 国产精品国产高清国产av| 国产成人精品无人区| 看黄色毛片网站| 听说在线观看完整版免费高清| 久久久久久久久中文| 又大又爽又粗| 日本黄大片高清| 色av中文字幕| 国产91精品成人一区二区三区| 中文字幕高清在线视频| 亚洲欧美一区二区三区黑人| 免费在线观看影片大全网站| 国产日本99.免费观看| 日韩欧美一区二区三区在线观看| 日日摸夜夜添夜夜添小说| 亚洲狠狠婷婷综合久久图片| 亚洲 欧美 日韩 在线 免费| 成年女人毛片免费观看观看9| 久久久国产欧美日韩av| 看免费av毛片| 欧美黑人巨大hd| 女同久久另类99精品国产91| 性欧美人与动物交配| 波多野结衣高清作品| 久久精品国产综合久久久| 高清毛片免费观看视频网站| 少妇人妻一区二区三区视频| 欧美日韩亚洲国产一区二区在线观看| 国产成年人精品一区二区| 少妇熟女aⅴ在线视频| 18禁黄网站禁片午夜丰满| 丁香欧美五月| 精品久久久久久久人妻蜜臀av| 啪啪无遮挡十八禁网站| 两性午夜刺激爽爽歪歪视频在线观看 | a级毛片在线看网站| 久久欧美精品欧美久久欧美| 又爽又黄无遮挡网站| 亚洲av成人av| 国语自产精品视频在线第100页| 欧美色视频一区免费| 国产三级在线视频| 50天的宝宝边吃奶边哭怎么回事| 丰满人妻一区二区三区视频av | 老鸭窝网址在线观看| 一区二区三区激情视频| 欧美极品一区二区三区四区| 国产91精品成人一区二区三区| 亚洲国产日韩欧美精品在线观看 | 日本免费一区二区三区高清不卡| 亚洲五月婷婷丁香| 2021天堂中文幕一二区在线观| 欧美中文日本在线观看视频| 国产av不卡久久| 91九色精品人成在线观看| 麻豆国产av国片精品| 日本免费一区二区三区高清不卡| 精品无人区乱码1区二区| 国产高清视频在线观看网站| 国产熟女xx| av福利片在线| 国产成人系列免费观看|