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

    Numerical investigation on the effect of depletion-induced stress reorientation on infill well hydraulic fracture propagation

    2022-03-30 13:52:32FngShouZhangLiuHuangLinYangEgorDontsovDingWiWngHongBoLiangZiRuiYinJiZhouTang
    Petroleum Science 2022年1期

    Fng-Shou Zhang ,Liu-K Huang ,c,* ,Lin Yang ,Egor Dontsov ,Ding-Wi Wng ,Hong-Bo Liang ,Zi-Rui Yin ,Ji-Zhou Tang

    a Key Laboratory of Geotechnical & Underground Engineering of Ministry of Education,Tongji University,Shanghai,200092,China

    b Department of Geotechnical Engineering,College of Civil Engineering,Tongji University,Shanghai,200092,China

    c State Key Laboratory of Geomechanics and Geotechnical Engineering,Institute of Rock and Soil Mechanics,Chinese Academy of Sciences,Wuhan,430071,China

    d ResFrac Corporation,Palo Alto,CA,94301,USA

    e PetroChina Research Institute of Petroleum Exploration and Development,Langfang,065007,China

    f State Key Laboratory of Marine Geology,Tongji University,Shanghai,200092,China

    g School of Ocean and Earth Science,Tongji University,Shanghai,200092,China

    Keywords:Infill well Depletion Hydraulic fracture Stress reorientation Asymmetry fracture

    ABSTRACT Depletion-induced stress change causes the redistribution of stress field in reservoirs,which can lead to the reorientation of principal stresses.Stress reorientation has a direct impact on fracture propagation of infill wells.To understand the effect of stress reorientation on the propagation of infill well's fractures,an integrated simulation workflow that combines the reservoir flow calculation and the infill well hydraulic fracturing modeling is adopted.The reservoir simulation is computed to examine the relationship between the extent of stress reversal region and reservoir properties.Then,the hydraulic fracturing model considering the altered stress field for production is built to characterize the stress evolution of secondary fracturing.Numerical simulations show that stress reorientation may occur due to the decreasing of the horizontal stresses in an elliptical region around the parent well.Also,the initial stress difference is the driving factor for stress reorientation.However,the bottom hole pressure,permeability and other properties connected with fluid flow control timing of the stress reorientation.The decrease of the horizontal stresses around the parent well lead to asymmetrical propagation of a hydraulic fracture of the infill well.The study provides insights on understanding the influence of stress reorientation to the infill well fracturing treatment and interference between parent and infill wells.

    1.Introduction

    Infill well fracturing technology is able to improve the reservoir's production efficiency.Figuring out the stress change caused by parent well's fracturing process and subsequent production is important for implementing infill well fracturing.The opening of the fracture and decrease of the pore pressure alter the stress field around the wellbore and fractures.Consequently,the direction of the maximum and minimum horizontal stressed is deviated from the original trend and the fracture propagation is influenced.Thus,it is vital to explore the mechanism behind the fracture propagation under the effect of stress reorientation.This topic is first studied by Mack and Elbel (1994) who find that along with the production,stresses within the depletion region decrease and the stress components in the direction parallel to the fracture decrease faster than the stress components in the direction perpendicular to the fracture.A reorientation theory is summarized based on this finding.This theory is improved by Siebrits et al.(1998) who bring up the concept of stress reversal region and analyze the influence of various reservoir properties by using a 3D simulation model.

    Fig.1.Concept of stress reversal region (modified from Siebrits et al.,1998).

    The concept of stress reversal region is illustrated in Fig.1.At the beginning,the direction of the maximum horizontal stress (σhmax)is parallel to the initial vertical fracture plane.During production,the stress parallel to the hydraulic fracture plane decreases faster than the stress in the direction perpendicular to it.Therefore,this alteration makes it possible to overcome the initial horizontal stress contrast between the minimum horizontal stress(σhmin)and σhmax.As a result,the stress reversal region develops around the fractures.In this case,refracturing at the same place could create a new fracture which initially propagates orthogonally to the initial fracture.Once it extends beyond the “isotropic point” at which the stresses in the two perpendicular directions are equal,the maximum principal stress switches back to its original direction,causing the secondary fracture to rotate back to be parallel to the initial fracture.The refracturing treatment makes use of this phenomenon so that the new fractures might intersect the region with more residual oil and higher fluid pressure (Elbel and Mack,1993;Zhang and Dontsov,2018).However,this phenomenon does not always occur in the field.The stress reorientation depends on many factors including the initial horizontal stress contrast,permeability,production rate and pore pressure drawdown (Li et al.,2006;Tan et al.,2020;Zeng et al.,2015;Zhai and Sharma,2007;Zhu et al.,2021).In order to better understand the hydraulic fractures’ deviation within the altered stress field,Berchenko and Detournay(1997) defined two dimensionless quantities:stress deviator and dimensionless time to minimize the number of independent parameters.Many simulation models are built to analyze the relationship among these factors and the extent of the stress reversal region,but few considered the coupled poroelastic mechanical effect on the stress reorientation (Aghighi et al.,2009;Chen and Teufel,2001;Lin et al.,2020;Roussel,2011).The opening of a propped fracture has been proved to be a source of the stress reorientation as well(Warpinski and Branagan,1989).Roussel and Sharma (2010) discussed the effect of poroelastic stress and stress induced by mechanical opening of hydraulic fractures together instead of independently and found that they are dependent on each other in a complex way,rather than being simply additive.

    One of the applications of the stress reorientation analysis is for infill well optimization.The key to a successful infill well fracturing is to make as large fracture surface area as possible and to penetrate less into the depleted zones.There are numerous previous works that focus on the fracture propagation within the perturbed stress field (Bruno and Nakagawa,1991;Gao et al.,2019;Li et al.,2017;Zhang et al.,2015).All these studies demonstrate that non-uniform pore pressure field does affect the initiation and propagation process of hydraulic fractures.Recently,Gupta et al.(2012)uncovered that reservoir depletion from parent wells can have detrimental impacts on fracture propagation of infill wells,and subsequently lowers the well productivity.In particular,in the case of low horizontal stress anisotropy,the intermediate stress direction is shown to reorient parallel to the infill well,potentially leading to longitudinal fractures.Not only production from offset wells impacts fracture propagation in the infill region,but stimulation of the infill well has also been shown to affect production in the offset wells,either positively or negatively(Ajani and Kelkar,2012;Gupta et al.,2012).In some cases,production in the offset wells may even be completely suppressed.Safari et al.(2017) proposed to take advantage of thein-situstress changes to prevent interference between the wells.Within the stress reversal region,direction of the maximum principal stress is deviated,which means that the infill well fractures may be deviated away from the preexisting wells.

    There is an extensive documentation of the pressure interference in shales aiming to achieve economic benefit of increased production from the infill well(Ajani and Kelkar,2012;Sahai et al.,2012).But in the case of infill drilling,hydraulic connection might not be the limiting factor.Whether poroelastic stress interference plays a bigger role than fluid flow interference depends largely on the value of thein-situstress contrast.In cases featuring large in situ horizontal anisotropy,stress reversal is unlikely to occur and fluid flow interference is then likely to dominate.Studying of the impact of pressure-induced stress changes on fracture propagation direction is by no means of recent interest (Berchenko and Detournay,1997).In tight gas plays,poroelastic stress interference has been previously investigated and applied to the refracturing of vertical wells (Roussel,2011;Roussel and Sharma,2010;Siebrits et al.,1998;Tan et al.,2021;Tang et al.,2019;Warpinski and Branagan,1989;Yi et al.,2020).

    Nevertheless,prediction of the behavior of hydraulic fractures within the perturbed stress field is a challenging task.The main objective of this study is therefore to investigate the propagation process of the fractures within the nonuniform stress field.However,it is common not to simulate the hydraulic fracturing process and production with the same model because these two processes have different nature and occur at different time scales.Therefore,the following workflow that combines the advantages of finite difference method and distinct element method is applied in this study.It is possible to obtain the stress field which has already been altered by the depletion in reservoir model by using finite difference method and model the hydraulic fracturing process by using the distinct element method with the input of the altered stress field.Thus,this study contains two parts and the first step is to examine the relationship between the extent of stress reversal region and reservoir properties.A set of dimensionless parameters is introduced to quantify the impact of mechanical properties,production parameters and geometry of the fractures to the stress reorientation.The sensitivity analysis of these variables is performed to illustrate the best practice of the infill well fracturing treatment.Then,the fracture propagation from a horizontal infill well and its sensitivity to stress reorientation is investigated.This work provides some insights on understanding the influence of stress reorientation to the infill well fracturing treatment and interference between parent and infill wells.

    2.Model setup

    2.1.Numerical scheme

    The finite difference method code FLAC3D is used to calculate the temporal and spatial evolution of the stress field in the vicinity of injection/production wells,which is able to simulate the rock deformation induced by pore pressure change(Itasca Consulting Group Inc,2011).Two main mechanical processes are included in this simulation:(1) the effective stress varies with the pore pressure,which affect the response of the rock,and (2) the fluid pressure reacts to volumetric deformations of the rock.The following equations considering Darcy's law,the fluid mass balance,mechanical constitutive laws and pore pressure are given to describe the hydro-mechanical coupling.

    The fluid transport in a porous material is described by Darcy's law:

    whereqiis the specific discharge vector,Pis pore pressure,ρfis fluid density,gkis gravity vector.i(j,k)take the values 1,2,3.kaijis apparent mobility coefficient,which is a function of fluid saturation.

    The fluid mass balance is formulated as:

    whereqvis the volumetric fluid source intensity.ζ is the variation of fluid content or variation of fluid volume per unit volume of porous material.

    The balance of momentum has the form:

    where σijis stress,is velocity,tis time.ρ=φρf+ρdis the bulk density,φ is porosity,and ρdis the density of the dry porous material.

    The variation of fluid volume is determined byPand ε (mechanical volumetric strains).The response equation for the pore fluid is expressed as:

    whereMis Biot modulus,and α is Biot coefficient.

    The elastic relations (constitutive law) that relate effective stresses to strains are formulated as (small strain):

    Where ˙εijis the strain rate,δijis the Kronecker delta,Δtis time step.KandGare bulk modulus and shear modulus of the porous material,respectively.The superscript0refers to the initial state.

    The relation between strain rate and velocity gradient(compatibility equation) is:

    The three-dimensional distinct element code 3DEC is employed in this study for the hydraulic fracturing modeling of an infill well,which is suitable for modeling a discontinuous medium that consists of blocks and discontinuities(Cundall,1988;Hart et al.,1988;Itasca Consulting Group Inc,2016).In this study,the rock mass is modeled as an assemblage of deformable rock blocks.Each individual block is discretized into an internal mesh consisting of tetrahedral zones that is used for finite difference calculations,while the joints represent interactions between the rock blocks through their boundaries.An explicit time stepping algorithm is employed for the mechanical calculation.At each time step,the constitutive equations are solved to calculate the stress in the whole domain.With the updated nodal force based on the updated stresses,the equations of motion are solved at all grid points to calculate velocities that are utilized to obtain the new strain rates.The integration of the equations of motion provides the new block positions and contact displacement increments.The joint forces are then updated based on the joint constitutive laws and applied to the blocks in the next time step.The fluid flow through the joints is modeled by the two-way hydro-mechanical coupled analysis in which the joint aperture and hydraulic conductivity depend on the mechanical deformation of the rock blocks,while the fluid pressure also affects the mechanical deformation of the rock blocks.The propagation of fluid-driven fractures is implemented by breaking the joints along the predefined joint plane in either tensile or shear mode(Huang et al.2019,2020;Yin et al.,2020;Zhang et al.,2019;Zhang and Mack,2017).

    2.2.Model geometry

    A reservoir production model is built to get the altered stress field by using FLAC3D.The model is implemented in three dimensions and contains a rectangular hydraulic fracture.The model configuration is illustrated in Fig.2 and the reservoir properties that are used for the base case are summarized in Table 1.

    Table 1 Geomechanical and reservoir properties for computation of the stress evolution due to depletion.

    The model consists of three layers including two bounding layers and one reservoir layer as the pay zone.The two bounding layers have the same heights and the layer of the reservoir has a variable height since the thickness of the pay zone may influence the extent of the stress reversal region (Siebrits et al.,1998).The reservoir properties are taken from the field case of Ordos Basin,China.The mechanical properties(Eb,νb)of the layers bounding the reservoir are different from the properties of the pay zone layer(Ep,νp).The flow is restricted to the reservoir layer by applying an extremely low permeability in the two bounding layers.The whole model is assumed to be elastic,homogeneous and isotropic except that the permeability is different in the horizontal and vertical directions.The boundary conditions for the model are illustrated in Fig.3 and can be summarized as:

    ●Fixed pore pressure at the face of the fractures

    ●Constant stress applied along the outer boundary of the domain:σzz=σv,σxx=σhmax,σyy=σhmin

    ●The outer boundary of the model is impermeable

    ●Fixed bottom hole pressure for production

    Considering the effect of boundary,the predefined hydraulic fractures are placed at a distance from the boundary equal to at least three times of the half fracture length in order to avoid the impact of boundary condition on the extent of the stress reversal region.Constantin-situstresses are applied at the outer boundary of the model.The bottom of the model is fixed.The hydraulic fracture is created by firstly splitting the grid points along the specified fracture surface and then by applying pressure to the two fracture surfaces.Consequently,a rectangular fracture volume is formed.Then,beam elements are used as proppant to make sure that the fractures remain open during depletion.

    Thein-situstresses in the hydraulic fracturing model are defined by importing the stress field calculated by reservoir model.The reservoir and surrounding rock matrix in the hydraulic fracturing model are assumed to be elastic,homogeneous,and isotropic.Fig.4 shows the 3D model constructed for the fracture propagation analysis.The properties used in the hydraulic fracturing model are the same as the ones used in the reservoir model.The hydraulic fracturing model has smaller dimensions compared with the model in the reservoir model in order to assure that the altered stress field from the reservoir model can cover the whole region.The hydraulic fracturing model contains one infill and one parent well.The two wells are parallel to each other and have a distance of 170 m between them.The fracturing fluid is assumed to be incompressible.

    Fig.2.The configuration of reservoir production model,a with a single hydraulic fracture,b with two hydraulic fractures.

    2.3.Dimensionless parameters

    The former research has proven that there are typically 4 types of dimensionless parameters affecting reservoir's stress reorientation:dimensionless time,dimensionless stress deviator,dimensionless fracture height ratio and dimensionless shear modulus ratio (Berchenko and Detournay,1997;Roussel,2011;Siebrits et al.,1998).

    It has been well known that the stress reorientation phenomenon is influenced by many parameters,such as the reservoir geomechanical properties,wellbore pressure,horizontal stress contrast,and the bounding layer properties.In order to minimize the number of relevant variables,Berchenko and Detournay(1997)defined two dimensionless parameters:τ,the dimensionless time;Π,the ratio of the initial deviation stress to the difference between initial reservoir pressure and well pressure.

    Dimensionless time:

    where τ is the dimensionless time,cis the fluid diffusivity,tis time,Lxfis the half length of the hydraulic fracture,kis permeability,μ is fluid viscosity,φ is porosity,Sis hydraulic storage,ν is Possison's ratio,andEis the Young's modulus.

    Fig.3.a Boundary conditions for the reservoir model (XY-plane) and b geometry of a vertical fracture in the reservoir model (XZ-plane).

    Fig.4.The configuration of the hydraulic fracturing model.

    Dimensionless stress deviator:

    where Π is the dimensionless stress deviator,P0is the initial reservoir pore pressure,andPwis the bottom hole pressure.S0is far-field stress deviator.σ*is characteristic stress.When bottom hole pressure is kept constant,P*equals difference between farfield pressure and bottom hole pressure.

    For 2D geometry,it is concluded that the production induced stress is scaled by the dimensionless stress deviator and depends only on the dimensionless time.But for the case of 3D geometry,two new dimensionless parameters are introduced.The numerical results by Siebrits et al.(1998) show that the layer bounding the reservoir affects stress reorientation.Thus,Roussel(2011)proposed two additional parameters to take into account the influence of the reservoir height and the mechanical characteristics of the bounding layers.

    Dimensionless fracture height ratio:

    where γ is the dimensionless fracture height ratio,hfis the height of hydraulic fracture,Lfis the length of hydraulic fracture.

    Dimensionless shear modulus ratio:

    where β is the dimensionless shear modulus ratio,Gbis the shear modulus of bounding layer,Gpis the shear modulus of pay zone.

    Initial maximum and minimum stresses both decrease during production due to depletion of the reservoir.Since the maximum stress decreases faster,at some distance from the fracture along they-axis the maximum and minimum stresses become equal.This corresponds to the isotropic point,whereby the extent of stress reversal region is defined by the distanceLyffrom the fracture to the isotropic point.Position of the isotropic point varies with the dimensionless time.And it is also clear that this distance is related to the fracture length.Therefore,it is better to use a dimensionless quantity to define the extent of the stress reversal region:

    and it is defined as a function of stress deviator Π,dimensionless fracture height ratio γ and dimensionless shear modulus ratio β.

    3.Numerical results of stress reorientation

    Fig.5 illustrates a typical stress reorientation pattern and provides a better perspective to the extent of the stress reversal region.Fig.5 shows a top view(z=0)of the simulation model,in which the well has been producing for τ=100 with the parameters specified in Table 1.The white dashed contour line agrees with the demonstration model of Siebrits et al.(1998) (Fig.1),but outside of the elliptical reversal region,there are still two wings,which represent a difference from the former research.The occurrence of stress reorientation comes from two aspects:mechanical stress reorientation and poroelastic stress reorientation,which will be discussed respectively in the following sections.

    3.1.Mechanical stress orientation

    Creation of a hydraulic fracture generates an additional stress that is mostly perpendicular to the fracture face.Fig.6 shows direction of the maximum horizontal principal stress after fracture opening and clearly demonstrates that open fractures can significantly alter the stress contrast and cause the direction of maximum horizontal principal stress to rotate 90°in the vicinity of the fracture.This phenomenon is called mechanical reorientation.For mechanical reorientation,the appearance of stress reversal region is caused by opening of the hydraulic fracture.The net pressure is directly applied on the face of the fracture and,as stated before,the hydraulic fracture is assumed as two planes.Therefore,the stress around the fracture increases and the stress parallel to the net pressure increase much more compared to the stress perpendicular to the direction of net pressure.As a result,the stress reversal region appears and because the change of stress is caused by the force applied on the fracture planar face,the stress reversal region is nearly rectangular.However,poroelastic reorientation is induced by the decrease of pore pressure and,after years’production,pore pressure field distributes in elliptical form.Thus,the stress reversal region is nearly elliptical,as shown in Fig.5.

    Fig.5.Contour of the maximum horizontal stress direction after τ=100 of production(90° represents the direction that is perpendicular to the fracture and 0° represents the direction that is parallel to the fracture.The white lines indicate the direction of maximum horizontal stress).

    The rotation angle of the maximum horizontal principal stress is 90°in this reversal zone,except the tip of the fracture.This is because the magnitude of the shear stress is small,as indicated in Fig.7.Consequently,directions of the horizontal principal stresses are along thex-axis andy-axis.Thus,it is found that within this reversal zone the maximum horizontal stress is along they-axis and the minimum horizontal principal is along thex-axis.On the other hand,the stresses in the direction parallel to the fracture(σxx)and perpendicular to the fracture(σyy)can represent the impact of mechanical reorientation.So,it is convenient to use the differences of σyyand σxxalong they-axis to quantify the degree of stress reversal.The generated stress contrast decreases as the distance from the fracture increases and the point at which the generated stress contrast equals to zero is the isotropic point.

    Extent of stress reorientation depends on the mechanical properties of reservoir and is more significant if the difference between the initial maximum and minimum horizontal principal stresses is small.At the beginning of depletion,Fig.8a shows that the additional force(Pnet)that is applied on the rock in the direction normal to the fracture face increases σyyand its influence on σxxis much smaller.Compressive stresses and strains have a negative sign in the model.So,if the increase of σyyis sufficiently large,the generated stress contrast will overcome the initial horizontal stress contrast to create a zone with abs(σyy)>abs(σxx).Therefore,the smaller the horizontal stress contrast is,the easier this alteration will be accomplished.In addition to the initial horizontal stress contrast,the reservoir's mechanical properties also affect stress reorientation.As shown in Fig.8b,a bigger stress contrast is generated in a reservoir with a lower Poisson's ratio.A low Poisson's ratio implies that the deformation in the direction parallel to the fracture is small compared to the deformation perpendicular to the fracture.Besides,Poisson's ratio also influences the production process which will be discussed in the next section.

    The extent of the stress reversal region is directly proportional to the fracture's dimensions.Several height values are tested in order to determine the relationship between the two variables.Fig.9 show that as the value ofhf/hr(hrequals the reservoir height)increases,the generated stress contrast increases too.However,the fracture heights are different in Fig.9a,but the height ratios are all equal to 1.0.Even if the height of the fracture increases for the constant height ratio,the generated stress contrast keeps the same value.

    3.2.Poroelastic stress reorientation

    Production from the reservoir generates nonuniform pore pressure in the formation and also alters the stress field due to the poroelastic effect.This phenomenon is called poroelastic stress reorientation.Our investigation aims to answer why stress reorientation happens and how this altering stress field evolves through time.

    First,the depletion induced pore pressure gradient is mainly responsible for this phenomenon (Fig.10).This nonuniform pore pressure field makes the stress field all decrease and the horizontal stresses in the direction perpendicular to the fracture and in the direction parallel to the fracture decrease in different speed.Based on Hooker's law and Biot's poroelasticity theory (Biot,1955),it's convenient to get the relationship between increment of stress and pore pressure,

    Fig.6.a Direction of the maximum horizontal principal stress around a single propped-open fracture.b Direction of the minimum horizontal principal stress around a single propped-open fracture.

    Fig.7.Contour of the shear stress around a single propped-open fracture.

    Fig.8.a The stress change induced by the opening of a hydraulic fracture,b effect of Poisson's ratio to the mechanical stress reorientation.

    Assuming that the vertical stress σzz=-σvdoes not change,it is easy to get the following relationship,

    After solving the above equations,we obtain,

    According to Eq.(14),the induced stress consists of two parts.One part is the change of pore pressure and the other part is deformation caused by the change of effective stress.The two components have opposite effect on the rock mass.The decrease of pore pressure is like unload to decrease the stress around the fracture.However,the deformation induced by fluid flow makes the rock mass more compacted,which means the effective stress arises.Especially,the deformation in the direction parallel to the fracture is different than that in the direction perpendicular to the fracture.Fig.11 shows the contour of the difference of strain incremental (Δεxx-Δεyy) around a single fracture at different time.There is an elliptical region that satisfies the condition(Δεxx>Δεyy).In other words,in this region the horizontal stresses in the direction parallel to the fracture decrease faster or increases slower than the horizontal stresses in the direction perpendicular to the fracture.Consequently,a stress reversal region occurs.

    Fig.9.a Generated stress contrast for different fracture heights and ratios of fracture height to reservoir height.b Generated stress contrast for different height ratios(hf/hr=0.5,0.75,1.0).

    Fig.10.Contour of the pore pressure for a single fracture producing for τ=100.

    Second,the expansion of this reversal region is time dependent.It can be seen from Fig.12a that the extent of the stress reversal region slightly decreases at early production time.This phenomenon can be explained by the action of opposing effects on the stress distribution between poroelastic and mechanical stress reorientation in the stress reversal region.Mechanical effect reorients stresses by applying more force to the fractures and,on the contrary,production induced stress reorientation is due to unloading(Fig.13).Therefore,the maximum principal stress is redirected when the production begins,which causes the stress reversal region to shrink (Fig.12b).The maximum stress will be fully reoriented within a fairly short time of production.Whereafter,the stress reversal region grows back with time.For large production times,poroelastic stress reorientation reaches a saturation plateau,so that there is a maximum value for the extent of stress reversal region during production.

    Except the case in Fig.12a,it has been studied that for different values of the stress deviator,fracture dimensions and dimensionless shear modulus ratio,the stress reversal region will have different sizes.The dimensionless time τ and ratio λ are used to illustrate the change of the stress reversal region extent over time.Fig.14 demonstrate the evolution of the stress reversal region versus dimensionless time for different values of stress deviator.It is apparent from the results that the smaller deviator is,the larger the stress reversal region will be.However,the weight of all factors in the stress deviator can be different.The sensitivity of the extent of stress reversal to other parameters is thus smaller compared to the stress contrast.Therefore,for a stress reorientation to occur,a small initial horizontal stress contrast is necessary.The smaller the initial stress contrast is,the easier for the production induced stress will be to overcome the difference between the initial minimum and the maximum horizontal principal stresses.

    As discussed previously,the stress deviator strongly affects the extent of the stress reversal region.Nevertheless,the factors included in the stress deviator have different impact on the curves of the dimensionless time versus extent of the stress reversal region.Fig.15 shows that the magnitude of the depletion induced horizontal stress contrast firstly increases and subsequently stabilizes and remains constant.As the initial horizontal stress contrast increases,the decrease of the curves can be slower or for large production times the curves can be almost flat.

    Pressure drawdown not only influences the extent of the stress reversal region,but also has effect on the timing of the occurrence of the maximum extent of the stress reversal region.One dimensional flow is known to be a good approximation for the flow in the middle of the fracture.In this case,solution for the diffusion equation that governs the pressure disturbance is,

    From the above equation,it is clear that different values of bottom hole pressure(or ΔP)can lead to different time instants,at which certain pressure disturbance magnitude reaches a pointx.In addition,different mobility values also affect the saturation time of the isotropic point.This explains why the extent of the reorientation region will expand quicker when the mobility is higher(Fig.16).

    One of the aims of infill well optimization is that the stress reversal region can protect the parent well from the interference caused by the infill well.Thus,as the stress reversal region becomes larger,it is more likely that the protection can be achieved.

    3.3.The case of multiple fractures

    The impact of fracture spacing to stress reorientation is also investigated for a multi-fracture scenario.To simplify the problem,the fractures are considered to be parallel to each other.Fig.17 show a single production well with two hydraulic fractures.Fig.17a demonstrates that the two fractures are spaced at a distance that is equal to one half of the fracture length while the spacing in Fig.17b equals to one eighth of the fracture length.Based on the maximum horizontal principal stress trajectories shown in Fig.17a,one can observe that the stress reversal region of these two fractures is separated.The gap between the two stress reversal regions makes it possible that the infill well fractures can propagate and intersect the producer well.Due to effect of stress shadow,the stress reversal region of each fractures reduces compared to the single fracture case.

    As shown in Fig.17b,because the spacing between the two fractures is reduced,a continuous stress reversal region forms,which can deviate away the approaching infill well fractures.Therefore,an appropriate spacing of the fractures can help producer to form a stress shield against the secondary fractures from those infill wells.

    4.Numerical results of infill well fracturing

    The distinct element method is used to simulate the fracture propagation process under the effect of stress redistribution in the infill well fracturing treatment.As discussed previously,the fracture propagation is not only affected by the variations of stress field,but it is also influenced by the fracturing fluid and the injection schedule.Since the main goal of this study is to explore the impact of stress reorientation on fracture propagation,the properties of the reservoir are assumed to be constant except that the stress field after production is different in each simulation.The perturbed stress field is calculated in the reservoir model and imported into the hydraulic fracturing model.

    Fig.18 illustrates different time instants of infill well facture propagation with the well spacing of 170 m(Fig.4),which includes the effect of stress reorientation.All the results show that when the fracture penetrates into the production zone,the width of the opening crack increases immediately and finally forms an asymmetric fracture which propagates away from the infill well along thex-axis.This is because the production causes an elliptical region with lower σyy,similar to the situation of fracture height growth into a low stress zone(Zhang and Dontsov,2018).The small value of the closure stress makes the wing of the hydraulic fracture located in the depleted zone to be longer than the one located in the undepleted zone.

    Fig.11.Contour of the difference of strain incremental (Δεxx-Δεyy) around a single fracture,after τ=0,50 and 100 of production.

    Fig.12.a The extent of the stress reversal region for λ versus τ with Π=0.0423,β=1 and γ=0.2,b reversal region's evolution from hydraulic fracture's opening to 160 h'production.

    Fig.13.Stress distribution curve for τ=0 and τ=100.

    5.Discussion

    One way to avoid the appearance of asymmetric fracture propagation is to adjust the inter-well spacing.Spacing the wells further away decreases pressure gradients and subsequently the risk of occurrence of asymmetric fractures in the infill well.As Fig.19 demonstrates,after the infill well is placed 310 m away from the parent well rather than 170 m,the infill well fracture propagates more symmetrically compared to the last case.Nevertheless,some pitfalls can be identified quickly.First,altered stress field's range is significantly impacted by time,as illustrated in Fig.14.So,the infill well drilling timing is connected to the parent well shutting time.Also,infill drilling is mainly considered in cases where well placement was initially dictated by leasing constraints.As it turns out,well spacing is not an easy parameter to adjust.But then other geomechanical considerations related to inter-well stress interference during stimulation would predominate as a result of the reduced well spacing(Roussel,2011;Roussel and Sharma,2010,2012;Siebrits et al.,1998).

    Fig.14.a The extent of the stress reversal region for λ versus τ with different Π,β=1 and γ=0.2.b The extent of the stress reversal region for λ versus τ with different bottom hole pressure and different Π,β=1 and γ=0.2.

    Fig.15.The stress contrast along y-axis at different times.

    Another way to decrease pressure gradients responsible for stress decrease or reversal in the infill region would be by shuttingin the producing offset wells for some period of time before infilling(Gupta et al.,2012) or pre-loading the parent well as a defensive approach (Whitfield et al.,2018).While we have already verified the high sensitivity of stress induced by pore pressure decrease to the pressure drawdown,it is still unclear how long of a shut-in or pre-loading period would be required to reverse the direction and magnitude of the principal horizontal stresses back to theirin-situstate and make the extent of the stress reversal region shrink.But considering the low permeability of shale reservoirs,it becomes uncertain that if the potential benefit from the stress redistribution would certainly prevail against the value lost from delaying production in the offset wells.Future work will be conducted to address these questions.

    The onset time of stress reversal has large implications on the optimization of the moment when infill operations should be carried on.We already saw that past this time,asymmetric propagation of fractures initiated from the infill well,and consequently the volume of reservoir contacted by the stimulation,will be negatively affected.Production impacts not only the direction of principal horizontal stresses,but also the value of the horizontal stress contrast.

    Fig.16.The extent of the stress reversal region for λ versus τ for different values of permeability.

    Previous studies showed that interactions between hydraulic and natural fractures are largely influenced by the stress regime in which the induced fracture propagates(Tang and Wu,2018;Zhang et al.,2017).If the difference between minimum and maximum horizontal stresses is small,an induced fracture will have a larger tendency to branch-out when intersecting natural fractures,thus forming more complex fracture networks.There is a window of opportunity for improved stimulation efficiency if fracturing the infill well around (ideally slightly before) the time when principal horizontal stresses switch directions.As a result,infill operations should not be considered as an after-thought but instead an integral part of the field development plan right from the start,taking into consideration both the risk of fracturing an infill well past the onset time of stress reversal,and the opportunity for improved fracture complexity if stimulation takes place before that critical time.

    6.Conclusions

    Fig.17.Contour of the stress reorientation angle after production(τ=100) around two propped-open fractures a with fracture spacing equal to half of the fracture length,and b with fracture spacing equal to one eighth of fracture length.

    Fig.18.Top:contour of the infill well fracture aperture at different time instants of injection;bottom:contour of YY-stress at different time instants of injection.

    Fig.19.Top:contour of the infill well fracture aperture at different time instants of injection,bottom:contour of YY-stress at different time instants of injection.The infill well spacing is increased to 310 m,compared with the model in Fig.18.

    This paper provides a case study to explore the time dependency of stress reorientation,as well as its behavior on problem parameters,such as threshold value of the horizontal stress contrast,the influence of stress deviator,dimensionless shear ratio,dimensionless fracture size and mechanical properties of the reservoir.The reservoir simulation is computed using the finite difference method,which considers both poroelastic effect and mechanical effect on the stress redistribution.The increase of pressure drawdown by decreasing bottom hole pressure and increasing mobility can expand the stress reversal region.The value of the initial horizontal stress contrast is the primary factor for stress reorientation.

    Fracture propagation during stimulation of the infill well is studied using the distinct element method.Results show that the decrease of the horizontal stresses around the parent well lead to asymmetrical propagation of a hydraulic fracture of the infill well.The width of the fracture in the vicinity of the parent well is also much larger compared to the width near the infill well.The asymmetrical propagation of the infill well hydraulic fracture confirms that it is necessary to consider the temporal and spatial effects when optimizing the completion of the infill well treatments.

    Acknowledgment

    The authors acknowledge the support provided by the Scientific Research and Technology Development Project of CNPC(Grant No.kt2017-19-01-1),the National Natural Science Foundation of China(Grant No.41772286,No.42077247 and No.42002271),and PetroChina Innovation Foundation(Grant No.2018D-5007-0202),and Project funded by China Postdoctoral Science Foundation (Grant No.2021T140514),and Open Research Fund of State Key Laboratory of Geomechanics and Geotechnical Engineering,Institute of Rock and Soil Mechanics,Chinese Academy of Sciences (Grant No.Z020009).The useful discussion with Branko Damjanac from Itasca Consulting Group is also greatly appreciated.

    亚洲内射少妇av| 亚洲av国产av综合av卡| 国产精品久久久久久久电影| 国产不卡av网站在线观看| 丁香六月天网| 欧美3d第一页| 午夜激情av网站| 中国美白少妇内射xxxbb| 午夜免费男女啪啪视频观看| 日韩电影二区| 亚洲欧美色中文字幕在线| 在线观看一区二区三区激情| 久久久久久久亚洲中文字幕| 国产高清三级在线| 一个人免费看片子| 亚洲精品美女久久av网站| freevideosex欧美| 成人影院久久| a 毛片基地| av在线老鸭窝| 成年美女黄网站色视频大全免费| 一本—道久久a久久精品蜜桃钙片| 国产免费视频播放在线视频| 亚洲精品456在线播放app| 亚洲国产看品久久| 欧美精品人与动牲交sv欧美| 精品久久国产蜜桃| 久久久久国产精品人妻一区二区| 成人综合一区亚洲| 交换朋友夫妻互换小说| 国产精品欧美亚洲77777| 亚洲精品国产av成人精品| 内地一区二区视频在线| 国产男女超爽视频在线观看| 日日撸夜夜添| av不卡在线播放| 日日啪夜夜爽| 国产精品久久久久久久电影| 午夜激情久久久久久久| 亚洲欧美中文字幕日韩二区| 午夜福利乱码中文字幕| 日韩av在线免费看完整版不卡| 我的女老师完整版在线观看| 热99国产精品久久久久久7| 极品少妇高潮喷水抽搐| 中文字幕人妻丝袜制服| 高清av免费在线| 又黄又粗又硬又大视频| 亚洲av在线观看美女高潮| 少妇高潮的动态图| www.色视频.com| 久久韩国三级中文字幕| av黄色大香蕉| 天堂中文最新版在线下载| 桃花免费在线播放| 亚洲综合精品二区| 精品一区二区三卡| 久久免费观看电影| 少妇的丰满在线观看| 日日爽夜夜爽网站| 亚洲精品一区蜜桃| 国产免费福利视频在线观看| 一级a做视频免费观看| 男人爽女人下面视频在线观看| 国产一区有黄有色的免费视频| 国产高清不卡午夜福利| 国产一区二区在线观看av| 国产免费视频播放在线视频| videos熟女内射| 免费高清在线观看日韩| 视频区图区小说| 免费大片18禁| 九九爱精品视频在线观看| 少妇高潮的动态图| 18禁裸乳无遮挡动漫免费视频| 啦啦啦视频在线资源免费观看| 最后的刺客免费高清国语| 国产在线一区二区三区精| 夫妻性生交免费视频一级片| 亚洲av成人精品一二三区| 大片免费播放器 马上看| 久久人人97超碰香蕉20202| 日本黄色日本黄色录像| 成人毛片60女人毛片免费| 满18在线观看网站| 欧美bdsm另类| 成人影院久久| 男男h啪啪无遮挡| 欧美成人午夜精品| 男女午夜视频在线观看 | freevideosex欧美| 精品久久久久久电影网| 国产一区二区在线观看日韩| 少妇人妻久久综合中文| 亚洲天堂av无毛| 日韩不卡一区二区三区视频在线| 国产一区亚洲一区在线观看| 久久ye,这里只有精品| 91精品国产国语对白视频| 一区二区三区乱码不卡18| 久久久久国产精品人妻一区二区| 男人舔女人的私密视频| 97在线人人人人妻| 人妻 亚洲 视频| 色94色欧美一区二区| 国产一区二区在线观看av| 国语对白做爰xxxⅹ性视频网站| 男女午夜视频在线观看 | a级毛色黄片| 伦精品一区二区三区| 久久久久精品久久久久真实原创| 2021少妇久久久久久久久久久| 精品一区二区三卡| 一区在线观看完整版| av卡一久久| 新久久久久国产一级毛片| 久久99一区二区三区| 亚洲欧洲日产国产| 国产欧美日韩一区二区三区在线| 亚洲av日韩在线播放| 有码 亚洲区| 麻豆精品久久久久久蜜桃| 久久久久久久久久成人| 国产成人免费观看mmmm| 亚洲精品第二区| 国产探花极品一区二区| 久久精品久久久久久噜噜老黄| 亚洲婷婷狠狠爱综合网| 国产午夜精品一二区理论片| 久久这里只有精品19| 国产亚洲精品久久久com| h视频一区二区三区| 汤姆久久久久久久影院中文字幕| 久久午夜福利片| 在线观看一区二区三区激情| 如日韩欧美国产精品一区二区三区| 国产av国产精品国产| av在线app专区| 久久99一区二区三区| 国产精品一区二区在线不卡| 18禁动态无遮挡网站| 九九在线视频观看精品| 黑人欧美特级aaaaaa片| 乱码一卡2卡4卡精品| 永久免费av网站大全| 中国国产av一级| 久久久久久人妻| 全区人妻精品视频| 人人澡人人妻人| 国产乱人偷精品视频| 日本黄大片高清| 久久久久人妻精品一区果冻| 在线观看人妻少妇| 大香蕉久久成人网| 国产精品欧美亚洲77777| 老司机亚洲免费影院| 热99久久久久精品小说推荐| 亚洲美女视频黄频| 51国产日韩欧美| 久久鲁丝午夜福利片| 色婷婷av一区二区三区视频| 欧美成人午夜精品| 中文字幕另类日韩欧美亚洲嫩草| 综合色丁香网| 精品亚洲成a人片在线观看| 亚洲精品aⅴ在线观看| 成人漫画全彩无遮挡| 18禁观看日本| 日韩一本色道免费dvd| av国产精品久久久久影院| 永久免费av网站大全| 久久午夜综合久久蜜桃| av线在线观看网站| 母亲3免费完整高清在线观看 | 国产免费一区二区三区四区乱码| 久久这里有精品视频免费| 日本色播在线视频| 曰老女人黄片| 欧美日本中文国产一区发布| 夫妻午夜视频| av在线app专区| 午夜日本视频在线| av不卡在线播放| 欧美激情 高清一区二区三区| 亚洲av成人精品一二三区| 精品人妻偷拍中文字幕| 亚洲色图综合在线观看| 秋霞伦理黄片| 一级毛片 在线播放| 久久精品夜色国产| 国产av精品麻豆| 国产精品欧美亚洲77777| 女人被躁到高潮嗷嗷叫费观| 国产永久视频网站| 成人无遮挡网站| 两个人免费观看高清视频| av福利片在线| 午夜激情av网站| 五月开心婷婷网| 亚洲内射少妇av| 纯流量卡能插随身wifi吗| 久久99一区二区三区| 天美传媒精品一区二区| 久久久精品免费免费高清| 婷婷色综合大香蕉| 99热全是精品| 午夜福利网站1000一区二区三区| 美女福利国产在线| 欧美+日韩+精品| 成人国语在线视频| 欧美亚洲日本最大视频资源| 成人国语在线视频| 丝袜在线中文字幕| 一级片免费观看大全| 亚洲欧美日韩卡通动漫| 国产精品国产av在线观看| 日韩精品有码人妻一区| 99热全是精品| 大话2 男鬼变身卡| 丝袜美足系列| 97超碰精品成人国产| av片东京热男人的天堂| 纯流量卡能插随身wifi吗| 亚洲精品aⅴ在线观看| 免费观看无遮挡的男女| 国产成人一区二区在线| 97在线视频观看| 韩国精品一区二区三区 | 狠狠精品人妻久久久久久综合| 男女免费视频国产| 亚洲婷婷狠狠爱综合网| 熟女av电影| 中国三级夫妇交换| 天天躁夜夜躁狠狠躁躁| 香蕉国产在线看| 99久久精品国产国产毛片| 人人妻人人爽人人添夜夜欢视频| 一级a做视频免费观看| 黑人巨大精品欧美一区二区蜜桃 | 亚洲国产成人一精品久久久| 久久国产精品大桥未久av| 亚洲成国产人片在线观看| 亚洲成人手机| 国产乱人偷精品视频| 狂野欧美激情性xxxx在线观看| 男女高潮啪啪啪动态图| 日韩欧美精品免费久久| 韩国精品一区二区三区 | 高清av免费在线| 在线观看国产h片| 韩国av在线不卡| 深夜精品福利| 久久这里只有精品19| av电影中文网址| 在线观看一区二区三区激情| 国产黄色免费在线视频| 亚洲国产最新在线播放| 人人妻人人爽人人添夜夜欢视频| av片东京热男人的天堂| 黄网站色视频无遮挡免费观看| 色94色欧美一区二区| 欧美精品国产亚洲| 欧美日韩亚洲高清精品| 日韩av在线免费看完整版不卡| 久久 成人 亚洲| 久久精品国产鲁丝片午夜精品| 亚洲一码二码三码区别大吗| 欧美少妇被猛烈插入视频| 久久久久久久国产电影| 边亲边吃奶的免费视频| 男男h啪啪无遮挡| 国产淫语在线视频| 成人国产麻豆网| a 毛片基地| 国产一区二区激情短视频 | 大香蕉久久成人网| 好男人视频免费观看在线| 久久这里有精品视频免费| 看非洲黑人一级黄片| 男人爽女人下面视频在线观看| 亚洲,欧美,日韩| 亚洲综合色惰| 香蕉精品网在线| 亚洲av日韩在线播放| 人人妻人人添人人爽欧美一区卜| 国产又爽黄色视频| 99热国产这里只有精品6| 久久久久精品性色| 大陆偷拍与自拍| 久久久久人妻精品一区果冻| 国产日韩欧美在线精品| 少妇人妻 视频| 国产黄色视频一区二区在线观看| 成人亚洲欧美一区二区av| 成人18禁高潮啪啪吃奶动态图| 亚洲三级黄色毛片| 亚洲av.av天堂| 欧美成人精品欧美一级黄| 男女边摸边吃奶| 欧美 亚洲 国产 日韩一| 国产成人精品福利久久| 精品少妇久久久久久888优播| 久久这里有精品视频免费| 中国三级夫妇交换| 精品亚洲乱码少妇综合久久| 日韩精品有码人妻一区| 欧美 日韩 精品 国产| 亚洲av中文av极速乱| 2021少妇久久久久久久久久久| 国产精品.久久久| 美女xxoo啪啪120秒动态图| 国产淫语在线视频| 免费久久久久久久精品成人欧美视频 | 下体分泌物呈黄色| 纵有疾风起免费观看全集完整版| 欧美人与性动交α欧美软件 | 国产片特级美女逼逼视频| 免费av不卡在线播放| 精品少妇久久久久久888优播| 蜜桃在线观看..| 亚洲欧洲日产国产| 日韩精品免费视频一区二区三区 | 女人被躁到高潮嗷嗷叫费观| 制服诱惑二区| 欧美3d第一页| 午夜av观看不卡| 精品国产国语对白av| 尾随美女入室| 亚洲av电影在线观看一区二区三区| 久久狼人影院| xxx大片免费视频| 十八禁高潮呻吟视频| 美女主播在线视频| 女人被躁到高潮嗷嗷叫费观| 久久精品国产自在天天线| 国产亚洲一区二区精品| 蜜臀久久99精品久久宅男| 国产精品一区www在线观看| 亚洲av日韩在线播放| 夜夜骑夜夜射夜夜干| 久久久国产精品麻豆| 精品少妇久久久久久888优播| 尾随美女入室| 视频中文字幕在线观看| 欧美日韩综合久久久久久| 国产成人a∨麻豆精品| 十八禁网站网址无遮挡| 亚洲 欧美一区二区三区| 日韩成人伦理影院| 两个人看的免费小视频| 免费在线观看完整版高清| 欧美另类一区| 另类亚洲欧美激情| 欧美精品av麻豆av| 蜜桃在线观看..| 国产片内射在线| 亚洲精品久久成人aⅴ小说| 国产精品人妻久久久影院| 97在线人人人人妻| 菩萨蛮人人尽说江南好唐韦庄| 韩国精品一区二区三区 | 日韩精品免费视频一区二区三区 | 男女午夜视频在线观看 | 人妻 亚洲 视频| 高清av免费在线| 久久精品熟女亚洲av麻豆精品| 草草在线视频免费看| 美女视频免费永久观看网站| 国产熟女午夜一区二区三区| 在线观看www视频免费| 国产精品国产三级专区第一集| 欧美xxxx性猛交bbbb| 1024视频免费在线观看| 激情视频va一区二区三区| 国产极品天堂在线| 激情视频va一区二区三区| 黄片无遮挡物在线观看| 中文字幕制服av| 婷婷色综合大香蕉| 熟妇人妻不卡中文字幕| 亚洲伊人久久精品综合| 男女午夜视频在线观看 | 欧美老熟妇乱子伦牲交| 两个人看的免费小视频| 免费大片黄手机在线观看| 久久青草综合色| 国产在线视频一区二区| 啦啦啦视频在线资源免费观看| 欧美日本中文国产一区发布| 亚洲欧美日韩卡通动漫| 人人妻人人澡人人爽人人夜夜| 少妇猛男粗大的猛烈进出视频| 岛国毛片在线播放| 丝袜美足系列| 香蕉丝袜av| 春色校园在线视频观看| 久久99一区二区三区| a 毛片基地| 1024视频免费在线观看| 国产精品久久久久久精品古装| 国产精品无大码| 婷婷成人精品国产| 久久这里只有精品19| 精品亚洲成国产av| 国产精品国产三级国产av玫瑰| 亚洲少妇的诱惑av| 久久久久视频综合| 精品一品国产午夜福利视频| 精品久久国产蜜桃| 欧美老熟妇乱子伦牲交| 啦啦啦视频在线资源免费观看| 亚洲成色77777| 精品国产国语对白av| 国产精品久久久久久久电影| 我要看黄色一级片免费的| 永久免费av网站大全| 黄色视频在线播放观看不卡| 久久免费观看电影| 性色av一级| 亚洲熟女精品中文字幕| 欧美成人午夜免费资源| 精品久久蜜臀av无| 免费看不卡的av| 久久热在线av| www.av在线官网国产| 在线观看免费高清a一片| 90打野战视频偷拍视频| 高清欧美精品videossex| 最近最新中文字幕大全免费视频 | 五月伊人婷婷丁香| 黑人巨大精品欧美一区二区蜜桃 | 18+在线观看网站| 亚洲国产精品一区三区| 看免费成人av毛片| 亚洲国产最新在线播放| 伊人亚洲综合成人网| 97精品久久久久久久久久精品| 欧美亚洲 丝袜 人妻 在线| 另类精品久久| 搡女人真爽免费视频火全软件| 97在线人人人人妻| 一级a做视频免费观看| 中国美白少妇内射xxxbb| 一区二区三区四区激情视频| 成人亚洲精品一区在线观看| 午夜老司机福利剧场| 午夜福利在线观看免费完整高清在| 国产片特级美女逼逼视频| 汤姆久久久久久久影院中文字幕| 欧美3d第一页| 美女国产视频在线观看| 男女无遮挡免费网站观看| 伦理电影免费视频| av有码第一页| 国产精品久久久久成人av| 成人毛片a级毛片在线播放| 免费黄色在线免费观看| 免费观看a级毛片全部| 丝瓜视频免费看黄片| 亚洲精品乱久久久久久| 80岁老熟妇乱子伦牲交| 亚洲欧美日韩卡通动漫| 制服丝袜香蕉在线| 一级黄片播放器| 日本色播在线视频| 亚洲少妇的诱惑av| 一级毛片电影观看| 在线天堂中文资源库| 精品国产国语对白av| 成人毛片60女人毛片免费| 亚洲精品第二区| 男人操女人黄网站| 新久久久久国产一级毛片| 午夜激情久久久久久久| 尾随美女入室| 久久精品国产综合久久久 | 蜜桃国产av成人99| www.色视频.com| 男人添女人高潮全过程视频| 丰满少妇做爰视频| 永久网站在线| 啦啦啦啦在线视频资源| 国产麻豆69| 夜夜爽夜夜爽视频| 亚洲情色 制服丝袜| 女人久久www免费人成看片| 日韩免费高清中文字幕av| 1024视频免费在线观看| 内地一区二区视频在线| 免费久久久久久久精品成人欧美视频 | 18在线观看网站| 国产精品久久久久久久久免| 日韩三级伦理在线观看| 宅男免费午夜| 丝袜美足系列| 91午夜精品亚洲一区二区三区| 中文字幕制服av| 777米奇影视久久| 欧美老熟妇乱子伦牲交| 午夜91福利影院| 欧美精品人与动牲交sv欧美| 另类精品久久| 日韩三级伦理在线观看| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 国产免费一区二区三区四区乱码| 久久久亚洲精品成人影院| 亚洲国产欧美日韩在线播放| 免费不卡的大黄色大毛片视频在线观看| 两个人免费观看高清视频| xxxhd国产人妻xxx| 免费观看在线日韩| 亚洲成人一二三区av| 亚洲综合色惰| 欧美老熟妇乱子伦牲交| av电影中文网址| 日产精品乱码卡一卡2卡三| 自线自在国产av| 色哟哟·www| 欧美日本中文国产一区发布| 久久久久网色| 26uuu在线亚洲综合色| 国产又爽黄色视频| 久久青草综合色| 最近最新中文字幕大全免费视频 | 精品久久蜜臀av无| 国产探花极品一区二区| 一本大道久久a久久精品| 只有这里有精品99| 99九九在线精品视频| 国产成人精品无人区| 妹子高潮喷水视频| 一本—道久久a久久精品蜜桃钙片| 欧美激情 高清一区二区三区| 咕卡用的链子| 91国产中文字幕| tube8黄色片| 观看美女的网站| 肉色欧美久久久久久久蜜桃| 黄网站色视频无遮挡免费观看| 桃花免费在线播放| 国产精品不卡视频一区二区| 午夜福利影视在线免费观看| 日本av免费视频播放| 久久久国产欧美日韩av| 人人妻人人添人人爽欧美一区卜| 亚洲国产精品一区二区三区在线| 女的被弄到高潮叫床怎么办| 啦啦啦在线观看免费高清www| 亚洲久久久国产精品| 99久国产av精品国产电影| 亚洲精品,欧美精品| av在线老鸭窝| 国产成人av激情在线播放| 国产精品久久久久久久电影| 婷婷色av中文字幕| 国产精品人妻久久久影院| 国产一区二区三区综合在线观看 | 中文欧美无线码| 免费黄频网站在线观看国产| 亚洲情色 制服丝袜| 亚洲欧美精品自产自拍| 在线亚洲精品国产二区图片欧美| 国产精品欧美亚洲77777| 女性生殖器流出的白浆| 日韩 亚洲 欧美在线| 亚洲美女视频黄频| 欧美变态另类bdsm刘玥| 王馨瑶露胸无遮挡在线观看| 最新的欧美精品一区二区| 免费黄网站久久成人精品| 国产精品欧美亚洲77777| 久久久久久久久久人人人人人人| 久久人人爽人人片av| 99热网站在线观看| 免费在线观看完整版高清| 亚洲av免费高清在线观看| 亚洲天堂av无毛| 夜夜爽夜夜爽视频| 国产精品久久久久久av不卡| 亚洲一区二区三区欧美精品| 午夜福利视频在线观看免费| 亚洲情色 制服丝袜| 在线看a的网站| 自线自在国产av| 久久影院123| 日本免费在线观看一区| 亚洲国产成人一精品久久久| 精品久久国产蜜桃| 99国产综合亚洲精品| 性色av一级| 国产男女内射视频| 久久精品人人爽人人爽视色| 一级黄片播放器| av又黄又爽大尺度在线免费看| 精品酒店卫生间| 少妇人妻精品综合一区二区| 亚洲av.av天堂| 国产一区二区三区综合在线观看 | 国产国语露脸激情在线看| av视频免费观看在线观看| 丰满少妇做爰视频| av免费观看日本| 热99久久久久精品小说推荐| 欧美 亚洲 国产 日韩一| 视频区图区小说| 亚洲精华国产精华液的使用体验| 99久国产av精品国产电影| 亚洲欧美成人综合另类久久久| 国产成人午夜福利电影在线观看| 欧美 亚洲 国产 日韩一| 视频区图区小说| 高清视频免费观看一区二区| 国产片特级美女逼逼视频| 亚洲综合色惰| 免费大片黄手机在线观看|