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

    Modeling of fiber bridging in fluid flow for well stimulation applications

    2020-06-22 06:04:48MehdiGhommemMustaphaAbbadGallyamAidagulovSteveDyerDominicBrady
    石油科學(xué)通報(bào) 2020年3期

    Mehdi Ghommem · Mustapha Abbad · Gallyam Aidagulov · Steve Dyer · Dominic Brady

    Abstract Accurate acid placement constitutes a major concern in matrix stimulation because the acid tends to penetrate the zones of least resistance while leaving the low-permeability regions of the formation untreated. Degradable materials ( fibers and solid particles) have recently shown a good capability as fluid diversion to overcome the issues related to matrix stimulation. Despite the success achieved in the recent acid stimulation jobs stemming from the use of some products that rely on fiber flocculation as the main diverting mechanism, it was observed that the volume of the base fluid and the loading of the particles are not optimized. The current industry lacks a scienti fic design guideline because the used methodology is based on experience or empirical studies in a particular area with a particular product. It is important then to understand the fundamentals of how acid diversion works in carbonates with different diverting mechanisms and diverters. Mathematical modeling and computer simulations are effective tools to develop this understanding and are efficiently applied to new product development, new applications of existing products or usage optimization. In this work, we develop a numerical model to study fiber dynamics in fluid flow. We employ a discrete element method in which the fibers are represented by multi-rigid-body systems of interconnected spheres. The discrete fiber model is coupled with a fluid flow solver to account for the inherent simultaneous interactions. The focus of the study is on the tendency for fibers to flocculate and bridge when interacting with suspending fluids and encountering restrictions that can be representative of fractures or wormholes in carbonates. The trends of the dynamic fiber behavior under various operating conditions including fiber loading, flow rate and fluid viscosity obtained from the numerical model show consistency with experimental observations. The present numerical investigation reveals that the bridging capability of the fiber– fluid system can be enhanced by increasing the fiber loading,selecting fibers with higher stiffness, reducing the injection flow rate, reducing the suspending fluid viscosity or increasing the attractive cohesive forces among fibers by using sticky fibers.

    Keywords Fiber bridging · Fiber flocculation · Modeling and numerical simulation · Discrete element method · Fiber– fluid coupling · Sensitivity analysis

    1 Introduction

    Accurate and efficient acid placement is a major concern in matrix stimulation of carbonate formations. During acidizing operations, the stimulation fluids will flow preferentially through the paths of least resistance, penetrating the zones with the highest porosity and permeability and then further increasing their conductivity while bypassing the lower-quality zones where stimulation is needed most.Moreover, in the presence of complex natural fractures and vugs in carbonates, the reservoir pressure and permeability may vary tremendously among zones. This makes achieving uniform stimulation and stable production levels a challenging task. To overcome these issues, diversion techniques have been developed with the aim of achieving effective placement of the stimulation fluids, creating a deep and uniform network with highly conductive flow paths (wormholes) between the reservoir and the wellbore,without causing reservoir damage (Asiri et al. 2013; Bukovac et al. 2012; Cohen et al. 2010; Garrouch and Jennings 2017; Ghommem et al. 2015; Kam et al. 2007; Qiu et al.2018; Siddiqui et al. 2006).

    Several alternatives are possible (and commercially available) to effectively place and divert stimulation fluids, and these can be grouped into two main categories: mechanical and chemical diversion (or the combination of both)(Alsaba et al. 2017; Alshubbar et al. 2018; Calcada et al.2015; Cohen et al. 2010; Davoodi et al. 2018; Droger et al.2014; Garagash et al. 2019; Gomaa et al. 2011; Kam et al.2007; Li et al. 2018; Malik et al. 2017; Sau et al. 2015;Solares et al. 2008; Thompson and Gdanski 1993; Wang et al. 2015a, b; Zhang et al. 2018). Mechanical diversion can be deployed for multistage stimulation completions. It is mainly performed by injecting ball sealers made of nylon,hard rubber or biodegradable materials such as collagen as part of a stimulation pumping sequence to plug and shut off the perforations and large fractures that are receiving most of the injected fluids (Alshubbar et al. 2018; Calcada et al. 2015; Cohen et al. 2010; Droger et al. 2014; Li et al.2017, 2018; Sau et al. 2015; Solares et al. 2008; Wang et al.2015a, b; Zhang et al. 2018). Alternatively, a straddle packer arrangement can be used to isolate the required treatment interval.

    Due to the completion and financial constraints, chemical diverters such as self-diverting acids (SDA) and viscoelastic diverting acid (VDA) (Schlumberger 2018b) are commonly used in carbonate acid stimulation. The use of these fluids has shown great potential for improving acid placement in carbonate reservoirs due to the rheological properties of the fluids that enable self-diversion to zones of lower injectivity(Asiri et al. 2013; Bukovac et al. 2012; Cohen et al. 2010;Garrouch and Jennings 2017). The SDA is a polymer-based in situ gelled acid. It consists of HCl mixed with a gelling agent and a pH-sensitive crosslinker. Initially, the SDA has a low viscosity. As the acid reaction with the formation takes place, the polymer crosslinks at pH = 1, which dramatically increases the viscosity. The resulting highly viscous acid slows down the fluid penetration rate and diverts the following acid stages to other intervals of the reservoir. At pH = 3.5, the gelled acid breaks and reverts to low viscosity. The narrow pH window limits the controllability of the SDA and therefore restricts its efficiency and its broad use.Moreover, polymer-based diverters leave residue after the gels break, potentially compromising the stimulation result by damaging the matrix or the surface of the wormholes,especially if the volume is improperly designed (Garrouch and Jennings 2017).

    To overcome the residue issue associated with polymerbased diverters, the viscoelastic surfactant-based diverting acid (VDA) was developed and introduced to the field(Garrouch and Jennings 2017). The viscosity increase in the VDA relies on the byproduct of the acid reaction with the formation rather than the fluid pH. The surfactant-based gel does not contain particles, which makes it nondamaging to the formation. However, if incorrectly applied or tested, the gel may not break in gas reservoirs or may form an emulsion with oil.

    In reservoirs that contain highly permeable streaks, such as natural fractures or vugs, the use of viscous gels may be not sufficient to provide efficient diversion, and fibers have been recently used as a diversion material. Degradable polymeric fibers are mixed with VDA to form a new generation of flow diversion agents (Schlumberger 2018a). The acid rapidly viscosi fies as it spends in the formation while the fibers form a porous medium that reduces or entirely blocks acid flow into thief zones (with high permeability)and redirects stimulation fluids to lower-permeability zones.Case histories of this new diversion fluid showed that the stimulation of highly heterogeneous carbonate reservoirs improved signi ficantly with the use of this novel acid-based fluid containing degradable fibers (Asiri et al. 2013; Bukovac et al. 2012; Cohen et al. 2010). For instance, as reported in a case study (Schlumberger 2011), despite performing five conventional acid stimulation treatments of a naturally fractured southern Mexico well (with a permeability contrast of more than 300:1) over a period of 3 years, a continuing decline in the production was observed. Furthermore, a production log showed that less than one-quarter of the interval was contributing to hydrocarbon production. Applying the degradable fiber diversion system provided a uniform zonal coverage, and the oil production stabilized at 254 m3/d, that is, 700% increase in the oil production.

    As presented above, there are several types of diverting agents for stimulating carbonates relying either on temporary plugging or viscosity control mechanisms (Alsaba et al.2017; Alshubbar et al. 2018; Asiri et al. 2013; Bukovac et al.2012; Calcada et al. 2015; Cohen et al. 2010; Davoodi et al.2018; Detournay 2016; Droger et al. 2014; Garrouch and Jennings 2017; Gomaa et al. 2011; Kam et al. 2007; Li et al.2018; Malik et al. 2017; Osiptsov 2017; Solares et al. 2008;Sau et al. 2015; Thompson and Gdanski 1993; Wang et al.2015a, b; Yu 2015; Zhang et al. 2018). Although these methods are often successful when applied during stimulation treatments, the industry still needs more reliable operational guidelines and a deeper understanding of the fundamental aspects of the diversion mechanism. The current empirical design methodology is based on experimental observations and mainly corrected by best practices. For this reason, an acid stimulation job may result in formation damage due to incorrect diverter volume, improper breaking and insufficient clean up energy. In these cases, the maximized poststimulation bene fit that cannot be ful filled is short lived or even not reached at all. Furthermore, when dealing with a fiber-based product, it is very hard for the operators to determine its merit and implement it with the highest con fidence.Therefore, an important step is to develop a qualitative and quantitative understanding of the interactions between the diverting agents and suspending fluids in the formation under different operating conditions. This understanding can be achieved by integrating advanced mathematical modeling and computer simulations with experimental work. This will help to recommend reliable diverting agents and operational guidelines for speci fic carbonate reservoir geology (lithology and heterogeneity).

    Fig. 1 Flow loop experiment. Left, general schematic. Right, photograph taken from Schlumberger Dhahran Carbonate Research Center (SDCR)laboratory

    The key mechanism behind the diverting capability of fibers is their ability to flocculate and bridge when interacting with suspending fluids. Since the introduction of fiber-based techniques for well intervention and stimulation applications, several laboratory-scale experimental studies(Alsaba et al. 2017; Alshubbar et al. 2018; Calcada et al.2015; Cohen et al. 2010; Li et al. 2018; Sau et al. 2015;Wang et al. 2015a) have been devoted to examining the fiber bridging phenomenon. The bridging experiments are conducted in pipe flow loops with slots of different shape and size mimicking fractures, as shown in Fig. 1. The bridging ability is tracked by monitoring the pressure signal and observing the fiber plug at the end of the experiment (Cohen et al. 2010; Sau et al. 2015). Li et al. (2018) proposed a novel fiber suspension system. They investigated the impact of the fiber length, fiber concentration, injection rate and permeability on the plugging efficiency. Their experimental study demonstrated the accessibility and good plugging effectiveness of the proposed fiber suspension system, and then they concluded that it can be efficiently applied in highpermeability reservoirs. Wang et al. (2015a) conducted a series of laboratory experiments using a triaxial stress system to investigate the diversion mechanisms of fiber-laden fluid systems for fracturing applications. They found that the dominant factors are mainly the fracture width, fiber concentration, heterogeneity and anisotropy of a reservoir,fluid viscosity and the injection rate. In a recent paper, Yang et al. (2019) performed experimental investigation to study the plugging mechanism of mixed fibers and particulates in a mimicked hydraulic fracture. To enable superior plugging performance, they found that fiber injection into the hydraulic fracture needs to be followed by injected particulates.Furthermore, the ratio between the fibers and particulates to be injected is suggested to be 1:1. Zhang et al. (2019) conducted a series of experiments on a 3D-printed rock model using a large-scale triaxial fracturing system to investigate the fiber-assisted diverting fracturing process. They identified the optimal fiber concentration from an economic perspective, fiber length and injection rate to enable temporary plugging and successful acid diversion.

    The present work focuses on the modeling aspect of fiber dynamics in fluid flow. The objective is to investigate the fundamental mechanisms of wellbore bridging by fibers of various characteristics when encountering openings (e.g.,natural fractures or wormholes) with different shapes and sizes. A transport model to predict fiber– fluid interactions in carbonates is developed and used to identify the key parameters in fluencing the bridging phenomenon and then affecting the performance of fluid placement.

    In this paper, we consider fluid diversion in a plug-andperf completion of a horizontal well with multistage acid fracturing. The idea is to inject fibers for temporary plugging previously formed fractures during acid fracturing operations. The goal is to control the uniform formation and growth of fractures over the entire well. To this end, we present a numerical investigation of the bridging performance of fibers in a mimicked fracture. Our study aims at integrating mathematical modeling and computer simulations with laboratory experiments to develop a fundamental understanding of fiber dynamics in fluid flow. Such understanding is intended to provide guidance on controlling efficiently the process of fiber bridging required to achieve accurate fluid placement in the wellbore.

    2 Mathematical modeling background

    2.1 DEM?based fiber model

    Due to the inherent complexity of fiber dynamics in fluid flow, particle-level simulations may be a useful complement to theoretical and experimental investigations. As such,many researchers have proposed fiber suspension models,in which the fibers are represented by multi-rigid-body systems of simple particles, for instance, interconnected spheres(as shown in Fig. 2) or ellipsoids. This approach is usually referred to as a discrete element method (DEM). One important feature of this method is the ability to model fibers with various equilibrium shapes, aspect ratio and stiffness. The geometry of the fiber can be then controlled by the size and the number of the particles (Elata and Berryman 1996; Guo et al. 2013; Switzer 2002; Switzer and Klingenberg 2003).

    Fig. 2 Discrete element-based representation of a flexible fiber

    The dynamics of fibers is modeled by tracking the motion of each particle (e.g., sphere) under internal interactions with the neighboring particles through bonds and external forces,such as hydrodynamic forces, contact forces with other fibers and gravitational forces. In summary, the dynamical process includes the following steps:

    · Detect contact among particles and their surroundings(e.g., walls).

    · Apply force–displacement law to compute the contact forces based on the relative motion between two entities at the contact and the contact model (as will be discussed in Sect. 2.2).

    · Apply the law of motion (Newton’s law) to each particle to update its velocity and position based on the resultant force and moment arising from the contact and any body forces acting on the particle. The governing equations obtained by applying the conservation of linear and angular momenta can be expressed as

    wherempandvpare the particle mass and velocity,Fsis the surface force,Fdis the force due to drag,Fpis the force due to pressure,Fbis the body force,Fgis the force due to gravity effect,Fconis the force due to contact developed between particles and their surroundings,Fcohis the attractive body force due to cohesion,Ipis the mass moment of inertia of the particle,ωpis the particle angular velocity andTcis the contact torque.

    2.2 Contact, bonding and cohesion

    The mechanisms that produce fiber aggregation are highly related to fiber– fiber and fiber–wall contacts. A contact model describes the physical behavior occurring at a contact(a small area where the physical touching and eventual overlap take place). Although a particle aggregation may exhibit complex dynamical behavior, this behavior may be achieved using relatively simple contact models. The contact model consists of up to three parts: (1) a contact-stiffness model,(2) a slip and separation model and (3) a bonding model.The contact-stiffness and the slip and separation models fully describe the physical behavior at all particle–wall and unbonded particle–particle contacts.

    A contact-stiffness model provides an elastic relation between the contact force and the relative displacement for unbonded particles. For instance, the damped linear spring as shown in Fig. 3 may be utilized. This model is composed of an elastic spring and a dashpot. The normal/tangential force (Fn/Ft) is generated by the repulsive elastic spring,and the damping force is applied by the dashpot, which dissipates energy during the contact. A more complex inelastic collision can be represented by the Hertz–Mindlin contact model (Elata and Berryman 1996).

    The slip model allows two entities in contact to slide relative to one another due to friction. They can also separate if a tensile force develops between them, and they are not bonded together. The slip condition occurs when the shear component of force reaches the maximum allowable shear contact force. This force is usually taken to be the minimum friction coefficient of the two contacting entities multiplied by the magnitude of the compressive normal component of the force.

    The bonding model makes the fibers elastically deformable. This model introduces attractive interparticle forces to the fiber. It uses the concept of a massless bar connecting a pair of bonded particles (see Fig. 4). The bar can transmit force and torque between neighboring particles and is subject to breakage under the maximum load. The bond radius controls the overall fiber flexibility. Connecting the spherical particles through larger bond radius yields fibers that are more rigid.

    Fibers in fluid flow may interact in two different ways:collisions and cohesive interactions. Collisions occur when the fibers are in physical contact and are usually described by contact-stiffness models discussed above. Cohesion is the property of particles (of the same material) to stick to each other due to mutual attraction. The cohesive interactions may be governed for instance by the square-well potential in which the attractive forces come into play after the separation distance between two particles reaches a speci fic distancedc.A typical interaction sequence is shown in Fig. 5. After the particles are close enough to each other, particles move in a straight-line trajectory, which leads to an inelastic collision.Because of the repulsive forces, particles travel in a straight path until the departing-cohesive interaction, which occurs at a separation distance ofdc. If the magnitude of these forces is not large enough, an internal re flection of the particles from the well occurs, which is referred to as a capture–cohesive interaction and yields to fiber aggregation. If the repulsive forces are signi ficant enough to overcome the cohesive well,an escaping cohesive interaction occurs and then the particles separate ( fiber dispersion). We note that for purpose of simplicity, the interaction in Fig. 5 is shown for only two particles and without incorporating other effects (e.g., wall, hydrodynamic). This physical aspect requires a careful investigation as being one of the key mechanisms of fiber flocculation and dispersion.

    Fig. 3 Schematic of a damped linear spring model

    Fig. 4 DEM interactions: particle bonding

    2.3 Fiber- fluid coupling

    The motion and interaction of fibers in a suspension/ floc/cake are strongly in fluenced by hydrodynamic forces exerted by the suspending fluid. As such, incorporating this aspect is important to analyze the bridging performance. The motion of an incompressible fluid phase in the presence of a secondary solid phase ( fibers) is governed by a modi fied set of Navier–Stokes equations given by

    Fig. 5 Fiber interaction sequence

    whereαfis the volume fraction occupied by the fluid,ρfis the fluid density,ufis the fluid velocity vector,τis the stress tensor of the fluid phase,pis the pressure, andCs?fdenotes the momentum exchange with the solid phase. We follow the Eulerian–Lagrangian approach to model the multiphase flow. As mentioned above, the Navier–Stokes equations are solved for the continuous phase (base fluid) while each formulated as described in Fig. 6. The geometrical and material properties of fibers along with the appropriate contact and cohesion models are first introduced. The geometry of the computational domain along with the fluid properties is speci fied as per the fracture and suspending fluid under investigation.individual particle within each fiber is tracked based on the Lagrangian approach. Such approach provides instantaneous information of the position and velocity of particles, and the forces being applied on them. It also has the capability to handle the flexibility of fibers, particle–particle collision(based on contact and cohesion models as described in the previous subsections), and the impact of walls. However, it may be computationally demanding for a large number of particles being tracked.

    Given the solvers for fluid motion (CFD-based) and fiber dynamics (DEM-based), and a coupling interface for the interaction between the two phases (e.g., through enforcing the conservation of momentum), a scheme for calculating the dynamics of fiber suspensions can be

    The coupling routines are detailed in the flowchart shown in Fig. 7. These include the following sequence of steps:

    · Clumped particles ( fibers) are injected into the fluid domain.

    · The position and orientation of each particle are obtained using the DEM-based approach and applying Newton’s law to each particle given the total force and torque that arise from the mechanical contact (developed between particles and their surrounding), bonding, cohesion and hydrodynamic effects.

    · For each particle, the corresponding mesh grid block is identi fied.

    · For each grid block, the volume fraction is determined.

    · The momentum exchange (source term in the modi fied set of the Navier–Stokes equations) between fluid and solid phases is calculated.

    · The flow and pressure fields are updated according to the current structure of the fiber suspension.

    · The hydrodynamic forces acting on the particles are computed and transferred to the DEM solver.

    · The same procedure is repeated. Fiber injection may be stopped when observing dominant fiber flocs and reaching the required fiber concentration. A careful selection of the timestep should be made to guarantee the convergence of the numerical solution (for both discrete fiber model and fluid flow) and avoid instability issues.

    Fig. 6 Modeling approach framework

    Fig. 7 Flowchart of the coupled computational model of the fiber– fluid system

    For the sake of simplicity and speeding up the numerical simulations, one can also employ the one-way coupling in which the impact of the fluid flow on fibers is accounted for but not the reverse. Other effects such as electrostatic, magnetic and heat transfer can be incorporated as well.

    Having a reliable coupled numerical fiber– fluid model,one may perform a sensitivity analysis of the tendency for fibers to flocculate when varying the fiber properties and the operating conditions. In particular, the size of the fiber flocs and their role in the bridging performance can be correlated with the characteristics of the fluid (viscosity, flow rate and regime) as well as with the fracture size. Such studies are limited (if not unfeasible) when relying only on experiments.

    Fig. 8 Fibers flocculated in suspensions. a Dry fibers. b Fibers dispersed in water-based fluid. c Entangled fibers forming a network structure.Pictures were taken from photographs courtesy of Schlumberger

    Fig. 9 Batch mixing of the degradable diverting acid system (Asiri et al. 2013)

    3 Fiber dynamics in fluid flow

    3.1 Fiber flocculation in sheared flow: mechanism of interest

    Flocculation describes the aggregation of fibers into clumps when a highly concentrated suspension of fibers is sheared.Figure 8 presents an example of flocs obtained from industrial fibers suspended in a water-based fluid. It shows clearly the capability of fibers to entangle, flocculate and form a network structure when interacting with fluid flow. Having such feature, fibers are viewed as an attractive option for acid diversion. During the acidizing jobs, operators use mixers (as shown in Fig. 9) to prepare the fiber-based diverting fluid (obviously by applying shear) prior to its injection to the well.

    Currently, there is a need to gain an insight into the fundamentals of how fluid diversion works in carbonates by different diverting mechanisms and diverters. As an attempt to address this, we develop a numerical model to examine fiber dynamics in fluid flow and investigate, in particular, the impact of different parameters on the bridging performance.

    Table 1 DEM-based fiber model properties

    3.2 Numerical model setup

    As a representation of a natural fracture, we consider a geometry consisting of a cylindrical pipe with a restriction at the end. This constitutes a representative geometry of lab-scale setup and imitating a natural fracture aperture to some extent. The variation in the pipe diameter before reaching the restriction is smooth (inclined boundary).The total length of the pipe is 5.5 in., and its diameter varies from 1 to 0.25 in. The dimensions of the simulated pipe with a restriction are designed to mimic the typical openings encountered in carbonates. At the inlet face(left of the pipe), we assume a uniform and constant inlet flow rate (Dirichlet boundary condition). A fixed pressure outlet (atmospheric pressure) is considered at the end of the pipe restriction (Neumann boundary condition). At the side walls, no penetration and no slip conditions are applied. The model geometry is built and meshed with the engineering software Star-CCM+ (Siemens 2017)in which all the computation and post-processing of the subsequent numerical results are performed. This software tool has demonstrated great capability to simulate multiphase flow with good accuracy in comparison with experiments (Marafaing et al. 2018; Tang et al. 2015). The system of incompressible fluid flow governing equations is numerically solved employing the finite volume method as a space discretization technique with an algebraic segregated solver. The geometry and computational mesh are shown in Fig. 10.As for the DEM-based fiber model, the particles are spheres with a diameter of 1.6 mm and a density of 1.4 g/L.More details on the particle properties and relevant model parameters are presented in Table 1.

    Fig. 11 Mechanical test to assess fiber rigidity. a Simply supported beam under central point loading and equivalent DEM-based discrete model.b Estimated fiber Young’s modulus as a function of bond radius multiplier

    3.3 Results and discussion

    3.3.1 Fiber calibration

    We evaluate the fiber Young’s modulus using a simply supported beam under a central point loading test (see Fig. 11a)while varying the radius of the bond that connects two neighboring particles and controls the fiber flexibility. Based on the Euler–Bernoulli beam theory, the maximum bending of the fiber can be expressed in terms of the fiber lengthL, the applied forceF, Young’s modulusEand the second inertia momentIas

    Table 2 Fiber– fluid system setup (numerical simulations shown in Fig. 12)

    whereris the particle radius. Therefore, by computing the central fiber bending, one could extract Young’s modulus.The variations of Young’s modulus with the bond radius multiplier (bond radius/particle radius) are shown in Fig. 11b. Clearly, increasing the thickness of the fibers of the bond connecting two neighboring particles within the fiber increases its stiffness. One could tune then the bond radius to obtain the fiber flexibility of interest. Note that Young’s modulus of the fibers simulated in the present study is in the order of few MPa.

    Fig. 12 Snapshots corresponding to successive timesteps during the flow of fiber flocs. Fibers are injected randomly from the inlet and undergo arti ficial attractive forces in the sphere placed near the inlet to form fiber clustering before reaching the restriction. These forces produce the flocculation behavior of the fibers when they get clustered to reproduce qualitatively the experimental observation as shown in Fig. 8

    3.3.2 Fiber clustering and bridging

    We perform a numerical experiment to study the fiber clustering and bridging in incompressible fluid flow. We first solve numerically for the fluid flow, check its convergence to steady state and then start injecting flexible fibers. Details on the relevant properties of the fibers and the suspending fluid are summarized in Table 2. The results are obtained while considering one-way coupling, that is, only the fluid flow affects the fiber motion and deformation. We present in Fig. 12 the set of snapshots corresponding to successive time steps during the fiber clustering and the flow of the formed fiber flocs along the pipe. As would be expected,the fluid flow is accelerated when reaching the restriction to ensure the conservation of mass. Color levels denote the size of the fiber clusters. Fibers are injected randomly from the inlet and undergo arti ficial attractive forces in the sphere placed near the inlet to form fiber flocs before reaching the restriction. Some of the fibers interact with the walls and slow down. The fibers of the first floc enter the restriction all together with a sliding motion between each other and lead to temporary and rapid bridging (t= 0.48 s).We observe a strong rearrangement of the fibers as they arrive at the restriction. The following fiber floc with bigger size (hitting the restriction at aboutt= 0.8 s) gets stuck at the restriction and leads to permanent bridging. The subsequently injected fibers continue to accumulate after bridging. In this case, the concentration of fibers is high enough to provoke many collisions between flowing fibers and prevent any rearrangement.

    3.3.3 Sensitivity analysis

    To characterize the fiber floc size, we introduce the mean floc diameter as a function of the fiber lengthL, the crowding numberNand the fiber aspect ratioc, as given by Switzer (2002):

    Hereφfis the fiber volume fraction. We define the floc–frac size ratio aswhereDis the restriction diameter. We note that the crowding parameter is considered in the present study to provide a better measure of fiber impact in comparison with the volume concentration of pure fibers (Lundell et al. 2011). Suspensions of flexible fibers in pipe flow are simulated for a variety of values of the system parameters. These include fiber Young’s modulus, mass flow rate and fluid viscosity. The effect of these parameters along with the floc–frac size ratio(proportional to fiber concentration) on the tendency for fibers to bridge is illustrated in Fig. 13. Flocs formed of rigid fibers tend to flocculate and bridge, whereas higher fiber concentration is required to achieve bridging when they are flexible. Being more rigid, they can resist more the pressure and drag forces at the pipe restriction, and then they can bridge more easily. This mechanism was observed experimentally; Soszynski and Kerekes (1988)reported that nylon fiber flocs dispersed when the fiber stiffness was reduced by heating above the glass transition temperature of nylon. Increasing the mass flow rate or fluid viscosity leads to higher drag and pressure forces acting on fibers and then this would provoke their rearrangements when reaching the restriction and prevent the bridging to take place. Thus, fiber flocs of larger size and higher concentration are needed to achieve bridging when dealing with more accelerated flow and/or viscous suspending fluids.

    Fig. 13 Effect of a fiber flexibility, b flow rate and c fluid viscosity on the fiber tendency to bridge

    Table 3 Key parameters affecting fiber tendency to bridge

    Table 3 summarizes the effects of different parameters on the tendency to bridge. The trends in the dynamical behavior of fibers under varying operating conditions are concluded based on numerical experiments. In addition to the effects discussed above, we examine the in fluence of the inter fiber cohesion (represented by a linear cohesion model in our numerical simulations, i.e., the application attractive forces proportional to the eventual particles overlap when colliding) on the fiber dynamics when interacting with the fluid flow. Activating the cohesion yields more compact flocs that promotes the tendency to bridge.

    Several models for bridging of particulates have been proposed in the literature (Garagash et al. 2019; Osiptsov 2017). In the pressure driven flow of a suspension through a fracture, the bridging is observed to occur when the ratio of the fracture width to the particle diameterw∕dis within the range of 2.5–3.0. A more recent study reported a modified version of the bridging criterion in a fracture which accounts for the particle volume fraction in the suspension flow (Osiptsov 2017). The critical fracture width is given by:

    where the bridging factorbis commonly taken equal to 2.5 andCpis the particle volume fraction. It was shown that the value of the bridging factor can vary signi ficantly based on the roughness of the walls (Garagash et al. 2019). The aforementioned bridging criteria currently used by several commercial simulation tools depend only on the particles size and concentration (Garagash et al. 2019; Osiptsov 2017). In the present study, we study the bridging mechanism of fibers in suspension flow. We found evidence of noticeable impact of the fiber stiffness, the base fluid viscosity and the inherent cohesion among fibers on the formation and stability of fiber flocs for temporary plugging.

    As a first step to verify the numerical predictions against experimental observations, based on flow loop tests conducted at Schlumberger Dhahran Carbonate Research Center, and establish a level of con fidence in the numerical transport model, we perform a set of simulations while reproducing most of the operating conditions ( fluid and fiber properties as reported in Table 4), and Fig. 14 shows the tendency for fibers to bridge at the restriction while varying the fiber loading and fluid flow rate. The shaded zone points out the range of values for the fiber concentration and flow rate that leads to fiber bridging. Clearly, increasing the flow rate lowers the tendency for fibers to bridge. The numerical results match closely with the experimentally observed fiber loading (15 g/L) needed to achieve bridging (represented by the dashed line in Fig. 14).

    Fig. 14 Bridging as a function of flow rate: experimental veri fication

    Based on the simulation results, recommendations are made for the hydraulic/acid fracturing operations in the field:(1) Stiffer fibers enable superior plugging performance at reduced fiber concentration/loading. (2) Adjusting the injection rate can be an alternative to optimize the fiber loading needed to achieve effective plugging from an economic perspective. (3) The viscosity of the base fluid needs to be accounted for when designing fiber-assisted diverting fracturing jobs given its signi ficant impact on the plugging mechanism.

    3.3.4 Fully coupled fiber– fluid system

    To obtain a more realistic representation of the fiber– fluid system, we perform a fully coupled simulation to examine the simultaneous effect of the fiber on the fluid flow and vice versa. To guarantee stability and convergence of the numerical solution, we reduce the time step to 10?4s.

    Table 4 Fiber– fluid system setup (numerical simulations shown in Fig. 14)

    We present in Fig. 15a, c the transient variations of the differential pressure (between the inlet and the outlet)and the corresponding snapshots of the flow of fiber flocs obtained from the numerical model. The first snapshot (see Fig. 15c) shows the flocculation of the fibers and their flow toward the restriction. At this level, the differential pressure is relatively low. As the fibers entangle and form a plug at the pipe restriction, the pressure rises. The pressure keeps increasing as the fibers accumulate and block the fluid to flow through the restriction until reaching a peak. The pressure buildup is the manifestation of the resistance of the bridge. Then, the pressure falls sharply to stabilize at a low value when the fibers exit entirely the pipe. The numerical results are in qualitative agreement with the experimental observations of the pressure trend when injecting fibers(Cohen et al. 2010) as shown in Fig. 15b, d. The flocculation and bridging process takes a few seconds. The experiments were conducted using the laboratory-scale equipment for testing leakoff behavior. As shown in Fig. 15d, this equipment includes a tube fitted with a piston, a pump and an ori fice. The acid and fibers are pushed through the piston while recording the transient pressure as the fiber-laden acid forms plugs and passes through the ori fice.

    Fig. 15 Transient variations of differential pressure while injecting fibers (upper frames) and flow of fiber suspensions (lower frames). Results are obtained from the numerical model and laboratory experiments for testing leakoff behavior. a Differential pressure versus time (numerical simulations). b Pressure signal during fiber-laden acid bridging experiment (Cohen et al. 2010). c Successive snapshots of the flow of a fiber suspension. Initially, the fibers are oriented randomly. They flip and rearrange as they flow toward the restriction and exit the pipe. d Experimental setup and example of fiber plugs formed in different geometry channels (Cohen et al. 2010)

    The computational time to simulate the dynamics of about 400 flexible fibers in fluid flow for 1-s physical time is approximately 12 h on 3.3 GHz processors. The simulation time is reduced to only 2 h when applying the one-way coupling and deactivating the impact of fibers on the fluid. As such, the fully coupled model is not suitable to conduct the parametric study which requires many simulations (hundreds if not thousands) to assess the bridging mechanism for different operating conditions. We have only shown few cases to demonstrate the possible tuning of the developed model to account for the two-way coupling. However, its use remains impractical for sensitivity analysis and optimization purposes.

    Fig. 16 Transient variations of differential pressure obtained from the fully coupled model (permanent fiber bridge)

    Next, we decrease the mass flow rate and investigate the bridging mechanism. Figure 16 shows the transient variations of the differential pressure resulting from fiber injection. The results are obtained from the fully coupled simulations for a physical time of one second. The sharp increase in the differential pressure is the manifestation of the formation of fiber bridge. The stabilization of the pressure level at higher value is indicative of the occurrence of permanent fiber bridge.

    Figure 17a shows the simulation time requirements of the case presented in Fig. 16. Figure 17b shows the corresponding simulation time when running the one-way coupled model. The results are presented in terms of elapsed time per time step. Clearly, the one-way coupled model is much less computationally demanding.

    4 Summary and conclusions

    In this paper, we developed a numerical model to simulate fiber flocculation and bridging in fluid flow. A discrete model of flexible fibers is coupled with a fluid flow solver to account for the inherent simultaneous interactions. This study aims at predicting the fiber dynamics when the fibers interact with suspending fluids and encounter restrictions that can be representative of fractures or wormholes in carbonates. Based on the present numerical analysis,we observed that the bridging capability of the fiber– fluid system could be enhanced by·

    Increasing the fiber loading

    · Increasing the fiber stiffness

    · Reducing the flow rate

    · Reducing the suspending fluid viscosity

    · Increasing the attractive cohesive forces among fibers(sticky fibers flocculate and bridge more easily)

    These trends in the dynamical fiber behavior (tendency for bridging under varying operating conditions) obtained from the numerical model show consistency with experimental observations. Furthermore, the fully coupled numerical simulation gives a pressure response (usually used to track the bridging) with features similar to that measured experimentally. This modeling capability can provide a baseline and guidance to improve the efficiency of existing fiber-based systems and to design new ones to fit the needs of new potential well intervention applications.

    Fig. 17 Simulation time requirements. a Two-way coupling. b One-way coupling

    Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons.org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

    在线免费观看的www视频| 黄色女人牲交| 成年女人看的毛片在线观看| 国产探花极品一区二区| 可以在线观看毛片的网站| 欧洲精品卡2卡3卡4卡5卡区| 一级av片app| 51国产日韩欧美| 禁无遮挡网站| 久久天躁狠狠躁夜夜2o2o| 一个人免费在线观看电影| 90打野战视频偷拍视频| 亚洲自偷自拍三级| 黄色女人牲交| 久久精品国产亚洲av香蕉五月| 精品久久久久久,| 此物有八面人人有两片| 我的老师免费观看完整版| 久久国产精品影院| 内射极品少妇av片p| av专区在线播放| 精品99又大又爽又粗少妇毛片 | 日本 欧美在线| 国产精品精品国产色婷婷| 亚洲av日韩精品久久久久久密| 亚州av有码| 国产精品,欧美在线| avwww免费| 18禁黄网站禁片免费观看直播| 国产av在哪里看| 97碰自拍视频| 亚洲精品在线观看二区| 好男人电影高清在线观看| 日韩欧美 国产精品| 动漫黄色视频在线观看| 欧美区成人在线视频| 免费av不卡在线播放| 成人高潮视频无遮挡免费网站| 国产在线男女| 免费av观看视频| 亚洲av中文字字幕乱码综合| 亚洲欧美日韩高清在线视频| 国产精品一区二区免费欧美| 日本三级黄在线观看| 国产精品永久免费网站| 免费黄网站久久成人精品 | 国产av在哪里看| 亚洲自拍偷在线| 久久人人精品亚洲av| 亚洲精品在线观看二区| 国产av在哪里看| 精品人妻1区二区| 人人妻,人人澡人人爽秒播| 精品人妻偷拍中文字幕| 日韩欧美 国产精品| 国内少妇人妻偷人精品xxx网站| 亚洲成人久久性| 1000部很黄的大片| 国内精品美女久久久久久| 欧美区成人在线视频| 日本黄色视频三级网站网址| 国产精品久久久久久久电影| 在线天堂最新版资源| 婷婷精品国产亚洲av| 婷婷精品国产亚洲av| 精品不卡国产一区二区三区| 麻豆成人av在线观看| 成年女人看的毛片在线观看| 麻豆久久精品国产亚洲av| 精品国产亚洲在线| 午夜两性在线视频| 午夜福利在线观看免费完整高清在 | 啪啪无遮挡十八禁网站| av视频在线观看入口| 女人十人毛片免费观看3o分钟| 亚洲aⅴ乱码一区二区在线播放| 91在线精品国自产拍蜜月| 亚洲av第一区精品v没综合| 两个人视频免费观看高清| 老鸭窝网址在线观看| 嫩草影院入口| 国产麻豆成人av免费视频| 精品欧美国产一区二区三| 久99久视频精品免费| 国产成人av教育| 很黄的视频免费| 亚洲欧美精品综合久久99| 成人亚洲精品av一区二区| 亚洲精品乱码久久久v下载方式| 丰满人妻熟妇乱又伦精品不卡| 久久国产精品人妻蜜桃| 啦啦啦韩国在线观看视频| 赤兔流量卡办理| 国产午夜精品论理片| 人妻制服诱惑在线中文字幕| 在线十欧美十亚洲十日本专区| 免费看光身美女| 国产精品1区2区在线观看.| 久久国产精品影院| 国产高清视频在线播放一区| 日韩人妻高清精品专区| 此物有八面人人有两片| 国产精品久久久久久亚洲av鲁大| 亚洲成人精品中文字幕电影| 久久精品久久久久久噜噜老黄 | 日本三级黄在线观看| 久久伊人香网站| 亚洲欧美日韩高清专用| 在线观看午夜福利视频| av在线天堂中文字幕| 两个人视频免费观看高清| 国模一区二区三区四区视频| 久久午夜亚洲精品久久| 狂野欧美白嫩少妇大欣赏| 99热这里只有是精品在线观看 | 午夜福利在线在线| eeuss影院久久| 丰满人妻一区二区三区视频av| АⅤ资源中文在线天堂| 色哟哟·www| 欧美日韩乱码在线| 97热精品久久久久久| 男人和女人高潮做爰伦理| 能在线免费观看的黄片| 97超级碰碰碰精品色视频在线观看| 男人狂女人下面高潮的视频| 直男gayav资源| 国内精品久久久久精免费| 韩国av一区二区三区四区| 久久精品国产亚洲av涩爱 | 欧美区成人在线视频| ponron亚洲| 久久国产乱子伦精品免费另类| av欧美777| 97超视频在线观看视频| 成年女人永久免费观看视频| 欧美性猛交╳xxx乱大交人| 亚洲国产精品成人综合色| 国产亚洲精品久久久com| 中文字幕精品亚洲无线码一区| 亚洲第一欧美日韩一区二区三区| 免费黄网站久久成人精品 | 国产精品久久久久久精品电影| 久久人人爽人人爽人人片va | 美女 人体艺术 gogo| 老司机午夜福利在线观看视频| 久久久久久久亚洲中文字幕 | 日本精品一区二区三区蜜桃| 我要看日韩黄色一级片| 99国产精品一区二区蜜桃av| 国产亚洲精品久久久com| 美女高潮的动态| 又黄又爽又刺激的免费视频.| 亚洲成av人片在线播放无| 国产精品一区二区三区四区免费观看 | eeuss影院久久| 亚洲18禁久久av| 国产精品女同一区二区软件 | 丰满人妻一区二区三区视频av| 麻豆国产97在线/欧美| 国产一区二区三区在线臀色熟女| 国产国拍精品亚洲av在线观看| 亚洲av成人av| 精品乱码久久久久久99久播| 97超视频在线观看视频| 最近视频中文字幕2019在线8| 亚洲av熟女| 日韩高清综合在线| 亚洲人成伊人成综合网2020| 99久久成人亚洲精品观看| 99久久精品一区二区三区| 韩国av一区二区三区四区| 麻豆成人av在线观看| 嫩草影视91久久| 热99在线观看视频| 变态另类丝袜制服| 欧美zozozo另类| 好男人在线观看高清免费视频| 国产精品亚洲av一区麻豆| 最近视频中文字幕2019在线8| 欧美zozozo另类| 亚洲人成电影免费在线| 蜜桃久久精品国产亚洲av| 乱码一卡2卡4卡精品| 一本综合久久免费| 亚洲 欧美 日韩 在线 免费| 亚洲美女黄片视频| 日韩欧美三级三区| 国内毛片毛片毛片毛片毛片| 最后的刺客免费高清国语| aaaaa片日本免费| 黄色丝袜av网址大全| 亚洲午夜理论影院| 精品国产亚洲在线| 日韩成人在线观看一区二区三区| 在线免费观看的www视频| 超碰av人人做人人爽久久| 禁无遮挡网站| 精品午夜福利视频在线观看一区| 亚洲av成人精品一区久久| 日韩国内少妇激情av| 夜夜看夜夜爽夜夜摸| 久久国产乱子伦精品免费另类| 88av欧美| 最近中文字幕高清免费大全6 | 免费av毛片视频| a级毛片a级免费在线| 色播亚洲综合网| 免费在线观看亚洲国产| 国产又黄又爽又无遮挡在线| 久久久久久久久中文| 亚洲国产高清在线一区二区三| 少妇人妻一区二区三区视频| 国产av不卡久久| 制服丝袜大香蕉在线| 麻豆成人午夜福利视频| 91麻豆av在线| 亚洲av中文字字幕乱码综合| 在线观看66精品国产| 免费看a级黄色片| 一二三四社区在线视频社区8| 一个人观看的视频www高清免费观看| 亚洲人与动物交配视频| 天堂影院成人在线观看| 亚洲一区二区三区不卡视频| 每晚都被弄得嗷嗷叫到高潮| 色噜噜av男人的天堂激情| 亚洲狠狠婷婷综合久久图片| 好男人在线观看高清免费视频| 亚洲人成网站在线播| 日日摸夜夜添夜夜添av毛片 | 精品不卡国产一区二区三区| 欧美在线黄色| 久久久久久国产a免费观看| 国产亚洲精品久久久com| 精品99又大又爽又粗少妇毛片 | 啦啦啦韩国在线观看视频| av天堂中文字幕网| 美女大奶头视频| 色哟哟·www| 日本一二三区视频观看| 色综合欧美亚洲国产小说| 波多野结衣巨乳人妻| 国产午夜福利久久久久久| 一区福利在线观看| 深爱激情五月婷婷| 欧美成狂野欧美在线观看| 级片在线观看| 欧美精品啪啪一区二区三区| 午夜福利在线在线| 女生性感内裤真人,穿戴方法视频| 亚洲五月天丁香| 国产单亲对白刺激| 小蜜桃在线观看免费完整版高清| 亚洲第一区二区三区不卡| 一级a爱片免费观看的视频| 欧美性猛交╳xxx乱大交人| 久久九九热精品免费| 国产久久久一区二区三区| 国产成+人综合+亚洲专区| 欧美+亚洲+日韩+国产| 久久久久久久久久黄片| 日韩高清综合在线| 欧美另类亚洲清纯唯美| 搡老熟女国产l中国老女人| 午夜福利高清视频| 在线免费观看的www视频| www.www免费av| 亚洲男人的天堂狠狠| 内射极品少妇av片p| 国产三级在线视频| 亚洲av日韩精品久久久久久密| 在现免费观看毛片| 最近最新中文字幕大全电影3| 欧美精品啪啪一区二区三区| 亚洲自拍偷在线| 美女xxoo啪啪120秒动态图 | 小说图片视频综合网站| 美女免费视频网站| 听说在线观看完整版免费高清| 51国产日韩欧美| 国产又黄又爽又无遮挡在线| 亚洲色图av天堂| 在线观看一区二区三区| 国产黄色小视频在线观看| 麻豆av噜噜一区二区三区| 国产探花在线观看一区二区| 亚洲成人免费电影在线观看| 国产高清激情床上av| 欧美在线一区亚洲| 99热这里只有是精品在线观看 | 黄片小视频在线播放| 一边摸一边抽搐一进一小说| 亚洲久久久久久中文字幕| 亚洲,欧美,日韩| 人人妻人人澡欧美一区二区| 国产蜜桃级精品一区二区三区| 天美传媒精品一区二区| 一级av片app| 天天一区二区日本电影三级| 日本与韩国留学比较| 亚洲无线在线观看| 51国产日韩欧美| 免费人成在线观看视频色| 午夜福利在线观看免费完整高清在 | 国产免费男女视频| 国产男靠女视频免费网站| 亚洲狠狠婷婷综合久久图片| 热99re8久久精品国产| 少妇的逼水好多| 亚洲 国产 在线| 亚洲狠狠婷婷综合久久图片| 免费在线观看亚洲国产| 9191精品国产免费久久| 婷婷精品国产亚洲av在线| 国产精品日韩av在线免费观看| 欧美日韩福利视频一区二区| 丁香六月欧美| 久久精品人妻少妇| 欧美一区二区国产精品久久精品| 非洲黑人性xxxx精品又粗又长| 天堂av国产一区二区熟女人妻| 日本免费一区二区三区高清不卡| 国产黄a三级三级三级人| 少妇人妻精品综合一区二区 | 国内精品久久久久久久电影| 国产激情偷乱视频一区二区| 麻豆久久精品国产亚洲av| 有码 亚洲区| 非洲黑人性xxxx精品又粗又长| 国产亚洲精品av在线| 男人舔奶头视频| 深爱激情五月婷婷| 男女视频在线观看网站免费| 桃色一区二区三区在线观看| 熟妇人妻久久中文字幕3abv| 亚洲成人中文字幕在线播放| 啪啪无遮挡十八禁网站| 亚洲av不卡在线观看| 黄色日韩在线| 宅男免费午夜| 男人的好看免费观看在线视频| av在线老鸭窝| 国产免费一级a男人的天堂| 国产精品日韩av在线免费观看| 国产美女午夜福利| 在线播放无遮挡| 深夜a级毛片| 欧美成人一区二区免费高清观看| 日韩欧美国产在线观看| 午夜福利在线观看免费完整高清在 | 夜夜看夜夜爽夜夜摸| 变态另类丝袜制服| 国产日本99.免费观看| 伊人久久精品亚洲午夜| 久久久成人免费电影| 亚洲五月婷婷丁香| 男女那种视频在线观看| 成人av一区二区三区在线看| 综合色av麻豆| 国产精品伦人一区二区| 网址你懂的国产日韩在线| 宅男免费午夜| av女优亚洲男人天堂| 日韩欧美在线乱码| 久久国产乱子免费精品| 国产精品久久视频播放| 中文在线观看免费www的网站| 欧美日韩黄片免| 亚洲av日韩精品久久久久久密| 97超级碰碰碰精品色视频在线观看| 99国产精品一区二区三区| 国内揄拍国产精品人妻在线| 嫁个100分男人电影在线观看| 午夜精品一区二区三区免费看| 亚洲 国产 在线| 我要搜黄色片| 日本在线视频免费播放| 亚洲第一区二区三区不卡| 亚洲国产欧美人成| 欧美丝袜亚洲另类 | 欧美乱色亚洲激情| 男女做爰动态图高潮gif福利片| 午夜精品在线福利| 中文字幕av在线有码专区| 五月伊人婷婷丁香| 日韩高清综合在线| av专区在线播放| 亚洲精品日韩av片在线观看| 亚洲黑人精品在线| 男女那种视频在线观看| 亚洲专区中文字幕在线| 超碰av人人做人人爽久久| 亚洲国产精品久久男人天堂| 国产一区二区三区视频了| 九色成人免费人妻av| 久久草成人影院| 国产老妇女一区| 精品久久久久久久久av| 亚洲第一欧美日韩一区二区三区| 在线观看免费视频日本深夜| 国产黄a三级三级三级人| 国产精品久久久久久久久免 | 午夜福利欧美成人| 欧美成狂野欧美在线观看| 91狼人影院| 久久久久久久久大av| 99精品久久久久人妻精品| 久久久久久久久久成人| 国产精品久久久久久精品电影| 国产人妻一区二区三区在| 人妻制服诱惑在线中文字幕| 美女黄网站色视频| 一级黄片播放器| 精品福利观看| 国产精品三级大全| 亚洲av不卡在线观看| 国产极品精品免费视频能看的| 国产一区二区三区在线臀色熟女| 欧美极品一区二区三区四区| 我的老师免费观看完整版| 亚洲av一区综合| 国产在视频线在精品| 日韩欧美精品v在线| 欧美xxxx黑人xx丫x性爽| 啦啦啦观看免费观看视频高清| 精品人妻视频免费看| 看黄色毛片网站| www日本黄色视频网| 精品欧美国产一区二区三| 国产熟女xx| 日韩欧美三级三区| 99热这里只有是精品在线观看 | 两个人视频免费观看高清| 九色国产91popny在线| 日本五十路高清| 俺也久久电影网| 中文资源天堂在线| 免费av毛片视频| 99热这里只有是精品在线观看 | 亚洲av熟女| 成人欧美大片| 一区福利在线观看| 日韩 亚洲 欧美在线| h日本视频在线播放| 色综合站精品国产| 日韩欧美国产在线观看| 免费看a级黄色片| 成年人黄色毛片网站| 国产伦一二天堂av在线观看| 88av欧美| 色综合欧美亚洲国产小说| 久久精品91蜜桃| 免费大片18禁| 怎么达到女性高潮| 亚洲成人中文字幕在线播放| 热99re8久久精品国产| 在线国产一区二区在线| 男女下面进入的视频免费午夜| 国产美女午夜福利| 看十八女毛片水多多多| 久久久久久久午夜电影| 精品久久久久久久久久久久久| 欧美+亚洲+日韩+国产| 99riav亚洲国产免费| 亚洲,欧美,日韩| 久久中文看片网| 特大巨黑吊av在线直播| 三级国产精品欧美在线观看| 免费人成视频x8x8入口观看| 午夜日韩欧美国产| 我的女老师完整版在线观看| 国产单亲对白刺激| 两个人的视频大全免费| 亚洲五月天丁香| 欧美乱色亚洲激情| 麻豆国产97在线/欧美| 欧美国产日韩亚洲一区| 亚洲自拍偷在线| av女优亚洲男人天堂| 久久国产精品影院| 欧美一区二区亚洲| 欧美最黄视频在线播放免费| 久久精品夜夜夜夜夜久久蜜豆| 一级黄片播放器| 国产极品精品免费视频能看的| 国产高清有码在线观看视频| 国产在线男女| 国产精品,欧美在线| 婷婷丁香在线五月| 大型黄色视频在线免费观看| 午夜免费激情av| 伊人久久精品亚洲午夜| 成年女人永久免费观看视频| 乱码一卡2卡4卡精品| 精品久久久久久久久亚洲 | 国产成人影院久久av| 亚洲欧美日韩东京热| 欧美一区二区国产精品久久精品| 又爽又黄a免费视频| 国内揄拍国产精品人妻在线| 美女 人体艺术 gogo| 精品一区二区三区人妻视频| 在线观看午夜福利视频| 免费看光身美女| 亚洲狠狠婷婷综合久久图片| 国产精品电影一区二区三区| 色综合亚洲欧美另类图片| 国产单亲对白刺激| 国产精品久久久久久亚洲av鲁大| 男人的好看免费观看在线视频| h日本视频在线播放| 一区二区三区四区激情视频 | 欧美绝顶高潮抽搐喷水| 国产一区二区三区在线臀色熟女| 99热这里只有是精品50| 香蕉av资源在线| 久久草成人影院| 亚洲精品一区av在线观看| 亚洲最大成人手机在线| 免费av不卡在线播放| 天美传媒精品一区二区| 可以在线观看毛片的网站| 国产精品免费一区二区三区在线| 午夜福利成人在线免费观看| 欧美日本亚洲视频在线播放| 久99久视频精品免费| 午夜福利在线观看免费完整高清在 | 五月伊人婷婷丁香| 国内精品久久久久久久电影| 99久久久亚洲精品蜜臀av| 99riav亚洲国产免费| 国内精品一区二区在线观看| 一个人免费在线观看电影| 欧美高清成人免费视频www| 亚洲精品色激情综合| 老女人水多毛片| 高清日韩中文字幕在线| 亚洲不卡免费看| 久久精品国产自在天天线| 麻豆av噜噜一区二区三区| 久久久久亚洲av毛片大全| 无人区码免费观看不卡| 亚洲成人久久性| 美女cb高潮喷水在线观看| 女同久久另类99精品国产91| 色播亚洲综合网| 又黄又爽又免费观看的视频| 日日干狠狠操夜夜爽| 国产精品乱码一区二三区的特点| or卡值多少钱| 在线国产一区二区在线| 国产亚洲欧美在线一区二区| 免费在线观看日本一区| 天堂动漫精品| 99热精品在线国产| 十八禁网站免费在线| 精品久久久久久久久久免费视频| 色综合亚洲欧美另类图片| 丰满人妻熟妇乱又伦精品不卡| 男女做爰动态图高潮gif福利片| 蜜桃久久精品国产亚洲av| 午夜a级毛片| 亚洲精品乱码久久久v下载方式| 男人舔女人下体高潮全视频| 亚洲久久久久久中文字幕| av欧美777| 如何舔出高潮| 午夜福利欧美成人| 精品日产1卡2卡| 搡老岳熟女国产| 国产精品亚洲av一区麻豆| 亚洲美女黄片视频| 国产单亲对白刺激| 婷婷色综合大香蕉| 99热这里只有精品一区| 一区福利在线观看| 别揉我奶头 嗯啊视频| 久久99热这里只有精品18| 99热这里只有精品一区| 午夜免费激情av| 夜夜看夜夜爽夜夜摸| 国产亚洲精品久久久com| 精品国产亚洲在线| 久久性视频一级片| 亚洲中文字幕一区二区三区有码在线看| www.www免费av| 全区人妻精品视频| 成年人黄色毛片网站| 亚洲成人久久性| 日韩 亚洲 欧美在线| 欧美又色又爽又黄视频| aaaaa片日本免费| 国产真实伦视频高清在线观看 | 精品人妻熟女av久视频| 99久久久亚洲精品蜜臀av| 亚洲欧美日韩高清专用| 每晚都被弄得嗷嗷叫到高潮| 成人性生交大片免费视频hd| 欧美+亚洲+日韩+国产| 在线看三级毛片| 深夜a级毛片| 欧美极品一区二区三区四区| 国产高清激情床上av| 中文字幕人妻熟人妻熟丝袜美| 精品乱码久久久久久99久播| 亚洲欧美日韩高清在线视频| 成年女人永久免费观看视频| 国产精品人妻久久久久久| 深爱激情五月婷婷| 亚洲欧美日韩东京热| 91久久精品国产一区二区成人| .国产精品久久| 亚洲欧美日韩高清专用| 神马国产精品三级电影在线观看|