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

    Peridynamic Modeling and Simulation of Fracture Process in Fiber-Reinforced Concrete

    2021-04-27 10:29:24ZhuangChenandXihuaChu

    Zhuang Chen and Xihua Chu,2,*

    1School of Civil Engineering,Wuhan University,Wuhan,430072,China

    2State Key Laboratory of Structural Analysis for Industrial Equipment,Dalian University of Technology,Dalian,116023,China

    ABSTRACT In this study,a peridynamic fiber-reinforced concrete model is developed based on the bond-based peridynamic model with rotation effect (BBPDR).The fibers are modelled by a semi-discrete method and distributed with random locations and angles in the concrete specimen,since the fiber content is low,and its scale is smaller than the concrete matrix.The interactions between fibers and concrete matrix are investigated by the improvement of the bond’s strength and stiffness.Also,the frictional effect between the fibers and the concrete matrix is considered,which is divided into static friction and slip friction.To validate the proposed model,several examples are simulated,including the tensile test and the three-point bending beam test.And the numerical results of the proposed model are compared with the experiments and other numerical models.The comparisons show that the proposed model is capable of simulating the fracture behavior of the fiber-reinforced concrete.Aft er adding the fibers,the tensile strength,bending strength,and toughness of the fiber-reinforced concrete specimens are improved.Besides,the fibers distribution has an impact on the crack path,especially in the three-point bending beam test.

    KEYWORDS Peridynamics;fiber-reinforced concrete;fracture mechanics;numerical simulation;three-point bending beam

    1 Introduction

    Fiber-reinforced concrete is a composite of concrete and fibers,which is wildly used in construction engineering [1–6].Conventional plain concrete has good compressive strength,but the tensile strength is low.And the brittle failure characteristic of the concrete structures limits its application.To improve the performance of the concrete,adding a small number of short fibers is a proper choice.There are several types of fibers applied in reinforcing plain concrete,including steel fibers,glass fibers,synthetic fibers,carbon fibers,and natural fibers.Fiber’s formation,length,volume fraction,and aspect ratio have a significant impact on the reinforcement effect.The fibers can restrain the microcrack generation,enable the crack bridging after the crack propagation,and relieve the stress concentration near the crack tip.After adding the fibers,many engineering properties such as strength,toughness,and durability of the concrete structure are improved.

    Numerous experimental studies are made to test the mechanical behavior of the fiberreinforced concrete [7–16].Although the experiment test can directly obtain the mechanical properties of the specimen,such as the tensile strength,Young’s modulus,and so on,it is expensive and time-consuming.Besides,the size effect of the experimental specimen affects the accuracy of predicting the mechanical properties of the engineering structures.The size of the engineering structure is much larger than the laboratory size.Considering these disadvantages,the researchers develop many numerical models to analyze the mechanical behavior of the fiber-reinforced concrete,such as the discrete element system [17],smeared crack model [18–20],cohesive zone model [21],micromorphic model [22],mixture theory [23],lattice model [24,25],finite element method (FEM) [26–32],and so on.Although these numerical methods can successfully predict the mechanical behavior of the fiber-reinforced concrete,there are still some challenges in dealing with the complex fracture propagation and the random fiber distribution.

    Peridynamics (PD),proposed by Silling [33],is a nonlocal version of the continuum mechanics.In peridynamics,the equilibrium equations are written in spatial integral form instead of spatial differential form.This special technology rids the difficulty of spatial derivatives near the discontinuity and makes the peridynamics more suitable for solving damage problems [34–44].In the peridynamics,the matter is discretized into a group of material points,and the interactions between points are through the bond.The original peridynamic model is called the bond-based peridynamics,where the bond acts like the spring [45].This spring-like interaction relationship leads to the limitation of Poisson’s ratio (μ=1/3 in plane stress,μ=1/4 in plane strain and 3D).To remedy this shortcoming,the researchers develop the state-based peridynamics [46,47]and modified bond-based peridynamics [48–51].In the state-based peridynamics,the interaction between material points depends on all the bonds connected to these two points.And the nonlocal deformation gradient is defined,which can correspond to the one in the continuum mechanics.However,the state-based peridynamics have some problems such as computing expensive and zero-energy mode.In the modified bond-based peridynamics,the more reasonable interaction relationships of the bond are developed,which consider the rotational degree and the tangential interaction of material points.

    In recent years,many numerical models based on the peridynamics for fiber-reinforced concrete are developed.Yaghoobi et al.[52] develop the fiber-reinforced concrete models by using the non-ordinary state-based peridynamics and the micropolar peridynamics [53],where the fibers are modelled by a semi-discrete method.Zhang et al.[54] simulate the crack propagation of a three-point bending beam test by using the bond-based peridynamics,where the randomly distributed PVA-fiber bonds describe the strengthening effect from fibers on the concrete matrix.Xu et al.[55] propose a fiber-reinforced concrete model to analyze the fracture behavior of prefabricated beams,where the different interactions and materials in the fiber-reinforced concrete are modelled by the different types of bonds and material points.These peridynamic models can successfully predict the crack path and the load-force curves of the fiber-reinforced concrete structures,but the interaction between fibers and the concrete matrix is oversimplified,and the frictional effect of the fiber is not considered.Besides,in paper [54],the bond-based peridynamics has the limitation of Poisson’s ratio,and the limited Poisson’s ratiosμ=1/3 andμ=1/4 are different from the Poisson’s ratio of concreteμ=0.2,which leads to the numerical errors in simulating the fiber-reinforced concrete.

    Based on the bond-based peridynamic model with rotation effect (BBPDR) [49],this work proposes a peridynamic model of fiber-reinforced concrete.The BBPDR considers the shear deformation of the bond and solves the limitation of Poisson’s ratio.Considering that the fiber volume fraction is low (0.25%–2% in common cases) and the scale of the fiber is smaller than the concrete matrix,the semi-discrete method is applied to model the fibers [52,53].In the proposed model,the fibers are not modelled independently,and the fiber reinforcement is implemented as the improvements of the bond’s strength and stiffness.Compare with the existing peridynamic fiber reinforced concrete models [52–55],the proposed model considers the strengthening effect and frictional effect of the fibers separately.The friction process starts when the local damage emerges.And the friction process is divided into static friction and slip friction.

    The main contribution of this work is developing a peridynamic fiber-reinforced concrete model based on BBPDR.Compared with the FEM fiber-reinforced concrete model,the proposed model is more suitable for fracture analysis of fiber-reinforced concrete.And the proposed model considers the frictional effect between the fibers and the concrete matrix,which is more reasonable for the fiber modeling than the existing peridynamic fiber-reinforced concrete models.The numerical results show that the proposed model can effectively simulate the crack propagation of fiber-reinforced concrete.And the influence of fiber content and the frictional effect is investigated through the tensile test and bending test.The numerical results indicate that the fiber content has a great impact on the strength and toughness,and the frictional coefficient mainly influences the toughness.

    The paper is organized as follows:Section 2 introduces the basic theories of peridynamics and BBPDR.In Section 3,the peridynamic fiber-reinforced concrete model is introduced,including the fiber modelling,strengthening mechanism,and frictional effect.Section 4 gives the numerical implementation of the proposed model and the ADR (adaptive dynamic relaxation) method.In Section 5,four numerical examples are presented to demonstrate the effectiveness of the proposed model,including the convergence study,the tensile test of a plate with a single fiber,the tensile test of the plate with circular notches,and the three-point bending beam test.In Section 6,the discussions and conclusions are made.

    2 Bond-Based Peridynamic Model with Rotation Effect(BBPDR)

    In peridynamic theory,the matter is discretized to a group of material points.These material points interact with each other through the bond within a finite distanceδ(δ=3.015·dxin common cases [56],dxis the discretization grid size),see Fig.1.The finite distanceδis called the horizon size.The interaction between material pointsiandjis called the pair-wise forcefij.The motion equations are written in the integral form in peridynamics:

    whereρis the mass density,xiis the location of pointi,(xi,t)is the acceleration of pointi,Hiis the peridynamic horizon,see Fig.1,fijis the pair-wise force of pointiacting on pointj,Vjis the volume of pointj,b(xi,t)is the body force density.

    In bond-based peridynamics,the pair-wise forcefijis dependent on the deformation of the bond between pointiandj.The relative displacement vectorηand initial relative position vectorξare defined as:

    Figure 1:The interaction between material points

    The pair-wise forcefijis linear with the bond’s stretchs:

    whereω(|ξ|)denotes the distance influence function and equals to 1 in this study,cis the micromodulus of the bond.The direction offijis parallel with the deformed bondξ+η,see Fig.2a.The determination of micro-moduluscis by comparing the distribution of strain energy density between the peridynamics and continuum mechanics.The value ofcis:

    whereEis Young’s modulus,μis the Poisson’s ratio,tis the thickness,δis the horizon size.

    Figure 2:The pair-wise force:(a) Bond-based peridynamics,(b) Bond-based peridynamic model with rotation effect (BBPDR)

    Because of the only one material parameter in bond-based peridynamics,there exists the limitation of Poisson’s ratio (μ=1/3 in plane stress,μ=1/4 in plane strain and 3D).To solve this problem,Silling et al.[46,47] develop the state-based peridynamic theories,including ordinary state-based peridynamics and non-ordinary state-based peridynamics.However,the state-based peridynamics are computing expensive and have the zero-energy mode.Considering the shortcomings of state-based peridynamics,the researchers develop the modified bond-based peridynamic models [48–51].These modified bond-based peridynamic models extend the relation between peridynamic pair-wise force and relative deformation,which considers the rotations of material points and the tangential peridynamic pair-wise force.

    The bond-based peridynamic model with rotation effect (BBPDR) is an extended model of the bond-based peridynamics [49].In this model,the bond can transfer the pair-wise forces in both axial and tangential directions,see Fig.2b.In isotropic elastic material,the number of material parameters in BBPDR is consistent with the one of continuum mechanics.And the numerical simulations of the paper [49] show that the limitation of Poisson’s ratio is solved.The axial pairwise forcefxand tangential pair-wise forcefyis linear with the bond’s stretches in axial and tangential direction,respectively:

    wherecanddare the micro-modulus in axial and tangential directions,uxanduyare the relative displacement components in axial and tangential directions.Similar to the method of determining the micro-modulus in bond-based peridynamics [45,56],canddare determined by comparing the distribution of strain energy density between the BBPDR and the continuum mechanics.In the plane stress problem,the values ofcanddare [49]:

    Comparing Eq.(8) with Eq.(6),it shows that the micro-modulus in bond-based peridynamics is identical with the axial micro-modulus in BBPDR.When the fixed Poisson’s ratioμ=1/3 is applied,the tangential micro-modulusdequals 0 and the BBPDR degenerates into the bondbased peridynamics.This comparison shows that the bond-based peridynamics is a special case of BBPDR.Although the theoretical and numerical analysis of the paper [49] show that the BBPDR can solve the limitation of Poisson’s ratio,considering that the tangential stiffnessdshould be greater than 0 to avoid the negative work of tangential pair-wise force,the selection of Poisson’s ratio expands to?1<μ<1/3.

    In peridynamics,the damage of material is described by bond breakage.The bond will break when its stretchsreaches its critical valuesc,and the material with such character is called PMB (Prototype microelastic brittle) material [45].Considering that concrete is a kind of brittle material,the assumption of PMB material is applied in this study.The scalar factorμ(t,ξ)is defined to show the bond’s breakage:

    wheret′is the time before the current timet.When the bond breaks,it cannot transmit interaction for material points anymore.scis obtained by comparing fracture energy calculated in peridynamics and the critical energy release rateGc.For simplicity,we only present the result in a plane stress problem [56]:

    whereGcis the critical energy release rate,cis the micro-modulus in bond-based peridynamics,tis the thickness,δis the horizon size.

    The damage value is defined as:

    whereφ(xi,t)is from 0 to 1,larger value means more broken bonds and more damage near the pointi.

    3 Fiber-Reinforced Concrete Model

    Considering that the fiber content is low,and the fiber scale is small compared with the concrete matrix,a semi-discrete model is applied to model the fibers,see Fig.3.Similar to the fiber modelling of paper [52,53],the fibers are not modelled independently,and the fiber reinforcement is implemented as the improvements of the bond’s strength and stiffness.Fig.3 shows that the bonds passed by the fibers are enhanced,including the bond’s stiffness and strength,see Fig.4.

    whereE1andE0are Young’s modulus of enhanced bond and normal bond respectively,ωf(|ξ|)is the influence function of a bond,which is chosen asωf (|ξ|)=dx/|ξ|in this study,the smaller the distance between the material points and the fibers,the stronger the strengthening effect of fiber on the material points,Efis the limit value ofE1,sc1andsc0are the critical stretch of enhanced bond and normal bond respectively,scfis the limit value ofsc1.

    In this study,the frictional effect between fibers and the concrete matrix is considered.The frictional process starts when the damage value of the material points near the fiber reaches limit valueφm:

    wherexiis the location of pointi,according to the paper [45],φmis chosen as 0.3,Gfis a group of material points,where the distance between the material points and the fiber is less thandf,see Fig.5.In this study,dfequals 0.8·dx.

    Figure 3:The bonds are enhanced by the fibers (δ=2dx)

    Figure 4:Fibers enhance the bond’s stiffness and strength

    Figure 5:The influence range of the fiber

    The frictional process is divided into two stages:Static friction and sliding friction,see Fig.6.During the static frictionthe frictional force densityfris linear with the relative displacement,which is calculated relative to the location when the friction process begins.During the sliding frictionthe frictional force densityfris a constant.The relation betweenfrandis:

    wherekfis the frictional stiffness,which is related to the interface properties between the fibers and concrete matrix,in this study,kfis calculated as:

    where dxis the discretization grid size,Cis a constant and equals 4000 (dx=1mm) and 2000(dx=0.5 mm),canddare the micro-modulus in BBPDR,see Eq.(8).is the sliding limit ofand equalssc·dxin this study.

    Figure 6:The relation between the frictional force density fr and relative displacement ur

    Besides,the frictional force density is limited by Coulomb’s law of friction.The frictional force density cannot be greater than the limit valueft.Whenfr≥ft,frequals toft.The limit value of frictional force density is calculated as:

    wherecfis the friction coefficient,nis the number of the enhanced bonds,fiis the pair-wise force of the enhanced bond,niis the tangential direction vector of the fiber,Vicalculates the compression between material point and fiber.The friction force density exists only when there is compression between the concrete and fibers.

    The fibers are distributed randomly by using the “rand” function in MATLAB,and random numbers determine the center locations and angles of the fibers.In the actual situation,the fiber volume fraction is used to represent the fiber content in concrete,which is defined asVf/V,whereVfis the total volume of fibers andVis the volume of the fiber-reinforced concrete.However,in the 2D plane stress problem,the fiber volume fraction calculated byVf/Vis relatively large because the thickness of the plate is small.Similar to the damage value of Eq.(11),the strengthening degree eiis used to represent the local fiber distribution near the region of pointi,which is defined as:

    whereRiis the number of the enhanced bonds,Niis the number of total bonds of the pointi.The large strengthening degree means the dense distribution of the fibers near the material point.

    4 Numerical Implementation

    In peridynamics,the material domain is discretized into an assemble of material points to solve the Eq.(1) numerically.And the integration symbol of Eq.(1) is replaced with the summation notation.In this study,the discretization scheme of the material points is the uniform grid.Eq.(1) is rewritten as:

    wherenis the number of the neighbors within the horizon,nijis the direction vector of the bond.

    There are two widely used time integration methods in peridynamic numerical models:Explicit method and ADR method.In this study,the ADR method is applied to obtain the quasi-static solutions of Eq.(19).The same method is used in the papers [52–54],and the detailed discussions of the ADR method are made in paper [57] and book [56].In the ADR method,the quasi-static solutions are obtained by adding the virtual damping force.The time step dt and the mass density are fictitious,dt equals 1.0 in this study.Eq.(19) is rewritten as:

    where M is the fictitious mass matrix and equals tomii·I,miiis the diagonal component of M and satisfied the inequality equation of Gerschgorin theorem [56],I is the identity matrix,U is the vector of displacements of all material points,C is the damping matrix and equals tocn·M,P is the vector of internal force density,b is the vector of body force density,X is the vector of the location of all material points.According to the book [56],miiis chosen as:

    where dt is the time step,tis the thickness of the plate,δis the horizon size,dxis the discretization grid size.

    The displacement vector U is obtained by applying the central time difference scheme:

    wherenrepresents thenth step,andcnis calculated as:

    Besides,there is a limit value ofcn.Whencncalculated by Eq.(25) is larger than 2.0,cnequals 1.9.Such limitation ofcnguarantees numerical stability,and the numerical error is small.

    5 Numerical Examples

    In this section,four examples are simulated to validate the proposed model and analyze the fracture behavior of fiber-reinforced concrete,including the convergence study,the tensile test of a plate with a single fiber,the tensile test of a plate with two semi-circle notches,and the threepoint bending beam test.In the first example,a static simulation of a square plate under tension is made to study the convergence of BBPDR.In the second example,a tensile test of the plate with a single fiber is simulated to validate the effectiveness of the proposed fiber-reinforced model.And in the following two examples,the specimens with randomly distributed fibers are tested.The results of the proposed model are contrasted with the numerical results of other models and the experimental results.

    5.1 Convergence Study

    A square plate under tension is tested to study the convergence of BBPDR,including the effects of discretization grid size dxand the horizon sizeδon the numerical results.The geometry and boundary condition of the plate is shown in Fig.7.The plate with dimension 100 mm×100 mm×1 mm is subjected to a displacement load at the top side,and its left and bottom side are fixed in x-direction and y-direction,respectively.Young’s modulusEis 210 GPa,and Poisson’s ratioμis 0.2.

    Figure 7:The geometry and boundary condition of the square plate under tension (unit:mm)

    Eq.(19) is solved by an implicit scheme,which is detailed discussed in paper [51].In the implicit scheme,the peridynamic integral equations are solved similar to the FEM,and the element stiffness matrix of a bond is obtained through the principle of virtual work.

    Fig.8 shows the x displacement distribution along the line y=50 mm (shown in Fig.7)with different discretization grid sizes (dx=1.0,2.0,and 5.0 mm),and the horizon sizeδequals 3.015·dx.Fig.9 shows the x displacement distribution along the line y=50 mm with different horizon sizes (δ=2.0,3.015,and 4.0 mm),and the discretization grid size dxis fixed as 1.0 mm.The numerical results are compared to the theoretical solution:

    wherexis the x coordinate,εyis the strain in y-direction and equals 0.001,μis the Poisson’s ratio.The relative error is calculated as:

    whereunare the numerical results.

    Fig.8 indicates that the numerical results converge to the theoretical solution as dxdecreases.Fig.9 shows thatδ=3.015·dxhas the highest accuracy,which is consistent with the suggestion of [56].Too large or too smallδleads to errors in calculating the element stiffness matrix of a bond.The largest error appears at the left boundary because the x displacement is close to 0 at the region close to the left boundary.In the following examples,δ=3.015·dxis applied in all the examples.In the plate with a single fiber,dx=1.0 mm is applied.And in the examples of crack propagation,dx=0.5 mm is applied for high precision in predicting fracture behavior.

    Figure 8:The x displacement along y=50 mm with different dx:(a) x displacement,(b) Relative error

    Figure 9:The x displacement along y=50 mm with different δ:(a) x displacement,(b) Relative error

    5.2 Tensile Test of a Plate with a Single Fiber

    In this section,the plate with a single fiber under tension is tested to validate the proposed fiber-reinforced model.The geometry and boundary condition of the plate is shown in Fig.10.The plate with dimension 100 mm×20 mm×1 mm is subjected to a displacement load at the right side.The displacement load increases linearly with the time steps.The 60 mm fiber is located at the center of the plate,and a weak region with a width of 10 mm is defined at the center of the plate,see Fig.10,where only the bonds connected to the points inside the weak region can break.The number of the material points is 2140,and the number of the total time step is 2000.Tab.1 lists the mechanical parameters of the plate,where theEfandscfare chosen greater than theE0andsc0to amplify the strengthening effect of fiber.A similar example with the same geometry condition is tested in the paper [52].

    Figure 10:The geometry and boundary condition of the plate with a single fiber (unit:mm)

    The force-displacement load curves of the plain concrete and the fiber reinforced concrete are shown in Fig.11.At the initial stage of loading (0–0.01 mm),the stress wave concentrates at the loading area,so the calculated reaction force remains zero.After the stress wave passes through the plate and reaches the left boundary,the calculated reaction force begins to rise.The tensile strength of the plate is increased by 11.5% and the ultimate strain is increased by 20%with the fiber reinforcement.Fig.12 shows that the damage initiates at the right boundary of the weak region and passes through the plate.After adding the fibers,the path of crack changes to the direction along with the fiber,see Fig.12b.Fig.13 shows the strengthening degree of the material points (Fig.13a) and the frictional force density (Fig.13b).The strengthening degree is larger when the distance between the fiber and the material point becomes closer.The positive direction of the frictional force is defined as the direction of the load,see Fig.13b.The frictional force densities are distributed on the material points close to the fiber,and the directions of frictional force densities are opposite on the two sides of the crack surface,see Fig.13b.Hence,the material points tend to move toward the crack surface,and the crack propagation near the fiber is prevented.

    Figure 11:The force-displacement load curves of the plates under tension

    Figure 12:The crack paths of the plates under tension:(a) Plain concrete,(b) Fiber-reinforced concrete

    Figure 13:The distribution of the strengthening degree (a) and the frictional force density along the fiber (b)

    The influence of the frictional coefficientcfis investigated and four different frictional coefficients are chosen,cf=0.1,0.3,0.7,0.9.The frictional coefficient describes the frictional characteristics of the fiber,which is related to the fiber formation and fiber type.Fig.14 shows the comparison of force-displacement load curves with different frictional coefficients.With the increase of frictional coefficient,the descending parts of the force-displacement load curves become slower,and the toughnesses of the plates become greater.However,the frictional coefficient has little effect on tensile strength.

    Figure 14:The force-displacement load curves of the plate with different frictional coefficients

    Besides,the proposed model is compared with the numerical model of fiber-reinforced concrete in paper [54].For simplicity,the numerical model of paper [54] is named as the random bond model in the following discussions.In the random bond model,the bond-based peridynamics is applied to calculate the pair-wise force between the material points,see Eq.(4),and the strengthening effect of the fiber is implemented by a randomly distributed PVA-fiber bond.This relation is similar to the fiber’s strengthening effect of the proposed model,see Fig.4.The stiffness and critical stretch of the PVA-fiber bond are higher than the concrete bond.And the fiber volume fraction is defined as:

    whereLfiberis the sum of lengths of all the PVA-fiber bonds,Lallis the sum of lengths of all the bonds.The fiber volume fraction determined by Eq.(28) is influenced by horizon sizeδ,andRvdecreases whenδrises.

    The tensile test setup is the same as Fig.10.Different from the fiber modelling of the proposed model,the fiber in the random bond model is represented by randomly distributed PVA-fiber bonds.Fig.15 shows the distribution of PVA-fiber bonds,and the volume fraction equals 20.26%,which is higher than the common cases (0.25%–2%).This example is only used to compare the proposed model and the random bond model,which is not of practical significance.The mechanical parameters of the random bond model are listed in Tab.2.The mechanical parameters of the proposed model are similar to Tab.1 and only the frictional coefficientcfis changed to 0.9.

    Figure 15:The distribution of PVA-fiber bonds in the random bond model

    Table 2:Mechanical parameters of the random bond model

    Fig.16 shows the force-displacement load curves of plain concrete and fiber-reinforced concrete.Since the random bond model is developed based on the bond-based peridynamics,the tensile test of plain concrete is also simulated by using the bond-based peridynamics.The two curves of the plain concrete are almost the same.And the descending part of the random bond model is characterized by brittle fracture.

    Figure 16:The force-displacement load curves of different models

    To further investigate the reinforcement effect of the random bond model,the plates with the different critical stretch of PVA-fiber bondscfare simulated.Fig.17 shows the force-displacement load curves with four differentscf(scf=1.0×10?3,1.2×10?3,1.5×10?3,1.8×10?3).The tensile strength increases asscfincreases,but the brittle fracture that occurs at the descending parts remains the same.Fig.18 shows the force-displacement load curves of the proposed model with four differentscf(scf=2.2×10?3,2.4×10?3,2.8×10?3,3.0×10?3).Because of the low content of the enhanced bond,see Figs.10 and 13a,in the proposed model,scfhas little effect on the tensile strength but slows down the descending part of the force-displacement load curve.The toughness of the plate becomes greater with higherscf.

    Figure 17:The force-displacement load curves of the random bond model with different scf

    Figure 18:The force-displacement load curves of the proposed model with different scf

    The results shown in Figs.16–18 indicate that the random bond model is capable of simulating the strengthening effect of fibers,but hard to simulate the improvement of toughness.Besides,the limitation of Poisson’s ratio in bond-based peridynamics may lead to numerical errors when simulating the fiber-reinforced concrete (μ=0.2).The proposed model is developed based on the BBPDR,which expands the selection of Poisson’s ratio (?1<μ<1/3).Both the strengthening effect and toughness improvement can be simulated by the proposed model,see Figs.11–14.And the semi-discrete fiber modelling shows the distribution of fibers more comprehensively.

    5.3 Tensile Test of a Plate with Two Semi-Circle Notches

    In paper [52],the tensile tests of a plate with two semi-circle notches are simulated to illustrate the effect of fiber volume fractions on tensile strength.The geometry and boundary condition of the plate is shown in Fig.19.The plate with two symmetry semi-circle notches is subjected to a displacement load.The mechanical parameters are shown in Tab.3.The selection of parameters and the geometry condition are referred to as paper [52].The number of the material points of the plate is 8153,and the number of the total time steps is 20000.

    Figure 19:The geometry and boundary condition of the plate with two semi-circle notches(unit:mm)

    Table 3:Mechanical parameters of the plate with two semi-circle notches

    As discussed in Section 3,the fiber volume fraction is replaced with the strengthening degree to characterize the fiber content.Fig.20 shows the distributions of fibers of the plate.Two different numbers of fibers (100 fibers and 200 fibers) are applied,and the fiber volume fractions are defined asαand 2α,respectively.The strengthening degree of the plate with a higher fiber volume fraction is distributed more intensively.

    Figure 20:The fiber distributions of the plate with a circular notch:(a) vf=α (100 fibers),(b) vf=2α (200 fibers)

    Fig.21 shows the force-displacement load curves of the plate with different fiber volume fractions.After adding the fibers,the tensile strengths of the plates are increased by 9.7% (vf=α)and 23.6% (vf=2α) relative to the plain concrete.Different from the brittle fracture of the plain concrete,the descending parts of the fiber reinforced concrete are characterized by the ductile fracture.Fig.22 shows the crack path of the plate with different fiber volume fractions,and the results are compared with the paper [52].The crack path is not a straight line after adding the fibers,but the damaged area is still in the center part of the plate.The fibers have a slight effect on the crack propagation but do not change its main direction.

    Figure 21:The force-displacement load curves of the plate with different fiber volume fraction

    Figure 22:The crack propagation of the plate with different fiber volume fraction:(a) Plain concrete,(b) vf=α,(c) vf=2α,(d) Fiber-reinforced concrete,vf=0.6% [52]

    Besides,the influence of fiber distribution on crack propagation is investigated.As reported in paper [9],the strengthening effect of fibers is maximized when the fibers are aligned perpendicular to the crack path.In this example,four different fiber distribution with vf=αare selected,named as (a),(b),(c),and (d).The strengthening degrees of four different fiber distributions are shown in Fig.23.Fig.24 shows the force-displacement load curves of the plate with different fiber distributions.The tensile strength of the plate is slightly influenced,and the fiber distribution mainly influences the descent part of the curves.The curves of distribution (b) and distribution(c) descend more slowly than the other two fiber distributions.This phenomenon results from the concentration of fibers in the center part of the plate.Figs.23 and 25 show that the strengthening degrees in the center part (|x|<2.5 mm) of distribution (b) and distribution (c) are higher than the other two fiber distributions.The plate with a higher strengthening degree in the center part shows greater toughness in the tensile test.Although the strengthening effect of fiber distribution and fiber agglomeration on concrete is complex,the strengthening degree defined by Eq.(18) can characterize the fiber reinforcement effect to a certain degree.Fig.26 shows the crack paths of different fiber distributions.The fiber distribution influences the local damage,but the damages still mainly concentrate in the center of the plate.

    Figure 23:The strengthening degree of the plate with four different fiber distribution,and the total strengthening degrees located at center region are shown in Fig.25:(a) 156.1,(b) 190.8,(c) 198.6,(d) 150.9

    Figure 24:The force-displacement load curves of the plate with four different fiber distributions

    Figure 25:The total strengthening degree of the material points located at the center region (|x|<2.5 mm)

    Figure 26:The crack paths of the plate with four different fiber distributions,and the total strengthening degrees located at center region are shown in Fig.25:(a) 156.1,(b) 190.8,(c) 198.6,(d) 150.9

    5.4 Three-Point Bending Beam Test

    This section discusses a three-point bending beam test.The geometry and boundary condition of the beam is shown in Fig.27.The final central deflection is 0.25 mm.The beam has a dimension of 160 mm×40 mm×0.5 mm with a notch of 18 mm.The number of the material points of the beam is 26019.The notch location is changed with c=0,20,40 mm.As the notch location changes,the crack pattern changes from the pure mode I to the mixed mode of Mode I and Mode II.The mechanical parameters are shown in Tab.4.The selection of the mechanical parameters is referred to as paper [54].Fig.28 shows the distribution of fibers and the number of fibers is 100.The number of total time steps is 40000.

    Figure 27:The geometry and boundary condition of the three-point bending beam test (unit:mm)

    Figure 28:The distribution of fibers of the three-point bending beam test

    Table 4:Mechanical parameters of the three-point bending beam test

    Fig.29b shows the crack path of the proposed model with c=0 mm.The crack pattern is the pure mode I,and the crack propagates along the centerline of the beam.Comparing with the results of plain concrete,see Fig.29a,the crack path has some perturbation because of the fibers.As discussed in the previous sections,the addition of fibers does not change the main direction of the crack path but influence the local damage distribution.The results match well with the experiments of paper [54],see Fig.29c.Fig.30 shows the deflection-force curves of plain concrete and fiber reinforced concrete.At the initial stage of the loading (0–0.01 mm),the calculated reaction force is zero or less than zero.Because the stress wave concentrates at the loading area,and the stress wave reflection occurs at the bottom boundary.After the stress wave reaches the bottom boundary (0.01–0.05 mm),the calculated reaction force begins to rise.The flexural strength of the beam is increased by 26.2% after adding fibers.When the curves approach the peak value (0.05–0.13 mm),the vibration of the curves becomes severe.And the same phenomenon appears during the descending part (0.13–0.25 mm).The crack influences the stress wave propagation,which leads to the oscillation phenomenon of the reaction force.Fig.31 shows the deflection-CMOD (Crack Mouth Opening Displacement) curves of plain concrete and fiber reinforced concrete.After the deflection reaches 0.08 mm,the CMOD of fiber reinforced concrete is lower than the one of plain concrete,which indicates that the fiber inhibits the crack opening process.

    Figure 29:The crack path of three-point bending beam (c=0 mm):(a) Plain concrete,(b) Fiberreinforced concrete,(c) Fiber-reinforced concrete [54]

    Figure 30:Deflection-force curves of three-point bending beam (c=0 mm)

    Figure 31:Deflection-CMOD curves of three-point bending beam (c=0 mm)

    When the notch location is off-center,the crack pattern is the mixed mode of Mode I and Mode II.Figs.32 and 33 show the crack propagation of the beams with c=20 mm and c=40 mm.In Fig.33,the crack path of the fiber reinforced concrete with c=40 mm is different from the experimental results.This difference results from the different local fiber distribution and the Poisson’s ratio.The experimental results of the paper [8] show that the oblique angles of cracks are different when the beam is prepared with different fiber volume fractions.And the numerical results of papers [53,54] show that the oblique angles of cracks are different from the experimental results when the fiber distributions are different.Fig.34 shows the crack path of fiber-reinforced concrete withμ=0.29,which is similar to the experimental results shown in Fig.33c.Poisson’s ratio influences the ratio of the axial stiffness and the tangential stiffness,see Eqs.(7) and (8).Considering that the effect of Poisson’s ratio on crack propagation is not the objective of this study,only the simulated results withμ=0.29 are shown.Figs.35 and 36 show the deflection-force curves of the three-point bending beams with c=20 mm and c=40 mm.The flexural strengths of the beams with different notch locations are increased by 52.2% (c=20 mm) and 25.6% (c=40 mm) after adding fibers.

    Figure 32:The crack path of three-point bending beam (c=20 mm):(a) Plain concrete,(b) Fiberreinforced concrete,(c) Fiber-reinforced concrete [54]

    Figure 33:The crack path of three-point bending beam (c=40 mm):(a) Plain concrete,(b) Fiberreinforced concrete,(c) Fiber-reinforced concrete [54]

    Figure 34:The crack path of three-point bending beam with Poisson’s ratio μ=0.29

    Figure 35:Deflection-force curves of three-point bending beam (c=20 mm)

    Figure 36:Deflection-force curves of three-point bending beam (c=40 mm)

    6 Discussion and Conclusion

    This paper proposes a peridynamic fiber reinforced concrete model based on the bondbased peridynamic model with rotation effect (BBPDR).The fiber is modelled in a semi-discrete way,and the fiber reinforcement is implemented by the improvements of bonds’strength and stiffness.In the proposed model,the friction interaction between fibers and the concrete matrix is considered,which is divided into static friction and slip friction.

    Numerical examples validate the proposed model’s effectiveness in modelling the fracture behavior of fiber-reinforced concrete.After adding the fibers,the tensile strength,flexural strength,and toughness of the concrete are improved.The fibers influence the path of crack propagation and the descending part of the load-deflection curve.The numerical results of the proposed model are compared with the experimental results and numerical results of other numerical models.And the comparisons show the effectiveness and validation of the proposed model in simulating the fracture behavior of fiber-reinforced concrete.

    In future work,the authors will investigate the corresponding relation between the numerical parameters and the physical properties of the fiber.And the optimization analysis of the fiber distribution and the hybrid fiber combination will be studied,which helps design the fiber-reinforced concrete with better performance.

    Acknowledgement:The authors thank the editors and the reviewers for their suggestions.

    Funding Statement:The authors are pleased to acknowledge the support by the National Natural Science Foundation of China through contract/Grant Nos.11772237,11472196 and 11172216,and to acknowledge the open funds of the State Key Laboratory of Structural Analysis for Industrial Equipment (Dalian University of Technology) through contract/Grant No.GZ19110.

    Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.

    又紧又爽又黄一区二区| 亚洲 欧美一区二区三区| 久久久精品国产亚洲av高清涩受| 国产精品偷伦视频观看了| 久久人妻福利社区极品人妻图片| 亚洲男人天堂网一区| 99精品欧美一区二区三区四区| 欧美黑人精品巨大| 亚洲自拍偷在线| 激情视频va一区二区三区| 欧美日韩瑟瑟在线播放| 成人永久免费在线观看视频| 午夜成年电影在线免费观看| 亚洲精华国产精华精| 在线免费观看的www视频| 精品人妻1区二区| 嫩草影视91久久| 一进一出抽搐gif免费好疼 | 国产亚洲精品第一综合不卡| 亚洲精品中文字幕在线视频| 在线观看免费视频日本深夜| 咕卡用的链子| 亚洲av片天天在线观看| 韩国av一区二区三区四区| 一级毛片高清免费大全| 免费一级毛片在线播放高清视频 | 91老司机精品| 最新在线观看一区二区三区| 亚洲欧美激情综合另类| 欧美人与性动交α欧美精品济南到| 国产熟女午夜一区二区三区| 成年人黄色毛片网站| 男人的好看免费观看在线视频 | 免费看十八禁软件| 久久久久久久久久久久大奶| 中亚洲国语对白在线视频| 亚洲七黄色美女视频| 狠狠狠狠99中文字幕| 亚洲九九香蕉| 日韩欧美免费精品| 免费高清在线观看日韩| 俄罗斯特黄特色一大片| 男女床上黄色一级片免费看| 日韩欧美一区二区三区在线观看| 精品国产美女av久久久久小说| 亚洲伊人色综图| 搡老岳熟女国产| 啦啦啦免费观看视频1| 久久久国产欧美日韩av| 日本五十路高清| 男人舔女人的私密视频| 日韩欧美国产一区二区入口| 午夜精品国产一区二区电影| 身体一侧抽搐| 亚洲一区二区三区色噜噜 | 天天影视国产精品| 亚洲成人久久性| 亚洲伊人色综图| 亚洲片人在线观看| 国产av一区二区精品久久| 日韩精品免费视频一区二区三区| 村上凉子中文字幕在线| 午夜影院日韩av| 国产xxxxx性猛交| 亚洲国产毛片av蜜桃av| 欧美在线黄色| 久久精品aⅴ一区二区三区四区| 免费av毛片视频| 国产蜜桃级精品一区二区三区| 亚洲片人在线观看| 黄片小视频在线播放| 国产91精品成人一区二区三区| 久久精品aⅴ一区二区三区四区| 国产精品国产av在线观看| 日本 av在线| 午夜a级毛片| 久热爱精品视频在线9| 女性被躁到高潮视频| 一级a爱视频在线免费观看| 成人特级黄色片久久久久久久| 999久久久精品免费观看国产| 精品国产超薄肉色丝袜足j| 欧美日韩黄片免| 国产欧美日韩一区二区三区在线| 韩国精品一区二区三区| 久久国产乱子伦精品免费另类| 丁香欧美五月| 男女下面进入的视频免费午夜 | 欧美激情高清一区二区三区| 18禁美女被吸乳视频| 亚洲av熟女| 国产黄色免费在线视频| 丝袜人妻中文字幕| av超薄肉色丝袜交足视频| 麻豆国产av国片精品| 亚洲人成电影免费在线| 操出白浆在线播放| 成人国语在线视频| 夫妻午夜视频| xxxhd国产人妻xxx| 亚洲一区二区三区不卡视频| 日韩免费高清中文字幕av| 精品少妇一区二区三区视频日本电影| 国产高清视频在线播放一区| 欧美日韩亚洲高清精品| 每晚都被弄得嗷嗷叫到高潮| 日韩成人在线观看一区二区三区| 成人三级黄色视频| 国产蜜桃级精品一区二区三区| 在线天堂中文资源库| 19禁男女啪啪无遮挡网站| 亚洲九九香蕉| a在线观看视频网站| 国产极品粉嫩免费观看在线| 精品无人区乱码1区二区| 老司机福利观看| 高潮久久久久久久久久久不卡| 一区二区三区国产精品乱码| 黑人操中国人逼视频| 在线十欧美十亚洲十日本专区| 亚洲精品中文字幕在线视频| 国产精品影院久久| 在线观看一区二区三区| 亚洲精品一二三| 午夜福利免费观看在线| 亚洲一码二码三码区别大吗| 国产精品偷伦视频观看了| 97超级碰碰碰精品色视频在线观看| 少妇的丰满在线观看| 免费久久久久久久精品成人欧美视频| 亚洲在线自拍视频| 精品国产国语对白av| 亚洲av五月六月丁香网| 精品国产超薄肉色丝袜足j| 男女午夜视频在线观看| 少妇被粗大的猛进出69影院| 99在线人妻在线中文字幕| 91av网站免费观看| 99riav亚洲国产免费| 国产极品粉嫩免费观看在线| 久久久国产一区二区| 最近最新中文字幕大全电影3 | 不卡一级毛片| 男人舔女人的私密视频| 国产一区在线观看成人免费| 18禁国产床啪视频网站| 黄片小视频在线播放| 国产精品av久久久久免费| 国产精品 国内视频| 久久久国产成人免费| 欧美中文日本在线观看视频| 久久人妻av系列| 国产成人精品久久二区二区免费| 老司机深夜福利视频在线观看| 女同久久另类99精品国产91| 99久久久亚洲精品蜜臀av| 午夜影院日韩av| 人人妻人人爽人人添夜夜欢视频| 国产精品自产拍在线观看55亚洲| 国产成人精品无人区| 一级作爱视频免费观看| 成人三级做爰电影| 亚洲精品久久午夜乱码| 亚洲va日本ⅴa欧美va伊人久久| 热re99久久精品国产66热6| 亚洲精品一区av在线观看| 久热这里只有精品99| 国产精品九九99| 久久亚洲真实| 热99re8久久精品国产| www.熟女人妻精品国产| 国产精品电影一区二区三区| 手机成人av网站| 亚洲专区中文字幕在线| av欧美777| 色综合站精品国产| 亚洲人成电影观看| 很黄的视频免费| 国内久久婷婷六月综合欲色啪| 看免费av毛片| 丝袜人妻中文字幕| 国产欧美日韩一区二区精品| 久久久久国产一级毛片高清牌| 99国产精品99久久久久| 老熟妇乱子伦视频在线观看| 99久久久亚洲精品蜜臀av| 国产又色又爽无遮挡免费看| 亚洲久久久国产精品| 久久久久久久精品吃奶| 一进一出抽搐动态| 成人18禁在线播放| 成人精品一区二区免费| 亚洲专区国产一区二区| 午夜精品国产一区二区电影| 国产一区在线观看成人免费| 国产一区二区激情短视频| av网站免费在线观看视频| 久99久视频精品免费| 亚洲精品国产一区二区精华液| 窝窝影院91人妻| 97人妻天天添夜夜摸| 女人高潮潮喷娇喘18禁视频| 精品乱码久久久久久99久播| 一夜夜www| 日本一区二区免费在线视频| 国产片内射在线| 美国免费a级毛片| 亚洲一区二区三区色噜噜 | 欧美日韩乱码在线| 99精品久久久久人妻精品| 99久久综合精品五月天人人| 精品高清国产在线一区| 女同久久另类99精品国产91| 亚洲国产精品999在线| 国产又爽黄色视频| 一区福利在线观看| 欧美av亚洲av综合av国产av| 一区二区日韩欧美中文字幕| 亚洲狠狠婷婷综合久久图片| 大香蕉久久成人网| 欧美日韩亚洲国产一区二区在线观看| 黑人操中国人逼视频| 国产成人欧美| 少妇被粗大的猛进出69影院| av在线天堂中文字幕 | 欧美大码av| 女性被躁到高潮视频| 亚洲成人久久性| 免费少妇av软件| 日韩三级视频一区二区三区| 免费看十八禁软件| 精品无人区乱码1区二区| 国产在线精品亚洲第一网站| 高清欧美精品videossex| 欧美+亚洲+日韩+国产| 精品国产一区二区三区四区第35| 天堂动漫精品| 精品福利永久在线观看| 两人在一起打扑克的视频| 日日爽夜夜爽网站| 精品国产超薄肉色丝袜足j| 亚洲片人在线观看| 亚洲七黄色美女视频| 色婷婷av一区二区三区视频| 50天的宝宝边吃奶边哭怎么回事| 国产伦一二天堂av在线观看| 不卡一级毛片| 在线观看一区二区三区| 嫩草影院精品99| 日韩高清综合在线| 久久这里只有精品19| 国产激情欧美一区二区| 一级a爱视频在线免费观看| 最好的美女福利视频网| 日本精品一区二区三区蜜桃| 啪啪无遮挡十八禁网站| 日本一区二区免费在线视频| 国产高清国产精品国产三级| 手机成人av网站| 免费少妇av软件| 国产有黄有色有爽视频| 日本撒尿小便嘘嘘汇集6| 国产激情欧美一区二区| 妹子高潮喷水视频| 久久精品国产亚洲av香蕉五月| 国产精品久久久久久人妻精品电影| 精品一区二区三区四区五区乱码| 两个人看的免费小视频| 国产精品99久久99久久久不卡| 美女午夜性视频免费| 亚洲七黄色美女视频| 精品日产1卡2卡| 老司机深夜福利视频在线观看| 亚洲成人精品中文字幕电影 | 中文字幕人妻熟女乱码| 欧美最黄视频在线播放免费 | 国产午夜精品久久久久久| 日韩高清综合在线| 国产人伦9x9x在线观看| av免费在线观看网站| 97超级碰碰碰精品色视频在线观看| 亚洲 欧美一区二区三区| 亚洲 国产 在线| 亚洲av日韩精品久久久久久密| 国产精品秋霞免费鲁丝片| 黑人操中国人逼视频| 婷婷精品国产亚洲av在线| 午夜福利在线免费观看网站| 日本vs欧美在线观看视频| 国产伦一二天堂av在线观看| 村上凉子中文字幕在线| 性色av乱码一区二区三区2| 精品欧美一区二区三区在线| 欧美性长视频在线观看| 久久国产精品男人的天堂亚洲| 日本免费a在线| 久9热在线精品视频| 久久久久精品国产欧美久久久| 欧美精品亚洲一区二区| 亚洲激情在线av| 男人操女人黄网站| 亚洲色图综合在线观看| 嫩草影视91久久| 国产有黄有色有爽视频| 久久草成人影院| 99国产精品一区二区三区| 亚洲国产毛片av蜜桃av| 久久精品国产亚洲av高清一级| 亚洲国产中文字幕在线视频| 黄色成人免费大全| 黄色女人牲交| 制服诱惑二区| 日日夜夜操网爽| 国产野战对白在线观看| 黄频高清免费视频| 欧美另类亚洲清纯唯美| 国产精品1区2区在线观看.| 精品乱码久久久久久99久播| 真人一进一出gif抽搐免费| 久久 成人 亚洲| 亚洲自拍偷在线| 啪啪无遮挡十八禁网站| 亚洲欧洲精品一区二区精品久久久| 侵犯人妻中文字幕一二三四区| 国产精品国产av在线观看| 国产精华一区二区三区| 真人做人爱边吃奶动态| 亚洲自拍偷在线| 午夜福利在线观看吧| 久久人妻av系列| 好男人电影高清在线观看| 亚洲第一青青草原| 视频在线观看一区二区三区| 人人澡人人妻人| 日韩欧美国产一区二区入口| 亚洲国产精品一区二区三区在线| 国产高清videossex| 亚洲精品中文字幕在线视频| av在线天堂中文字幕 | 国产免费男女视频| 一a级毛片在线观看| 视频在线观看一区二区三区| 中文字幕av电影在线播放| 欧美+亚洲+日韩+国产| 狠狠狠狠99中文字幕| 黑人欧美特级aaaaaa片| 老汉色∧v一级毛片| 99精国产麻豆久久婷婷| 校园春色视频在线观看| 国产亚洲av高清不卡| 黄色视频不卡| 亚洲人成伊人成综合网2020| a级毛片黄视频| 中文字幕人妻丝袜一区二区| 老司机午夜十八禁免费视频| 欧美精品啪啪一区二区三区| 欧美老熟妇乱子伦牲交| 中文字幕人妻熟女乱码| 两性午夜刺激爽爽歪歪视频在线观看 | 亚洲一区中文字幕在线| 久久精品国产清高在天天线| 午夜精品久久久久久毛片777| 亚洲中文日韩欧美视频| 国产一区在线观看成人免费| 91国产中文字幕| netflix在线观看网站| 久久精品aⅴ一区二区三区四区| 午夜精品在线福利| 热re99久久国产66热| 在线观看一区二区三区| 露出奶头的视频| 免费看十八禁软件| 国产激情欧美一区二区| 亚洲国产精品999在线| 一边摸一边抽搐一进一小说| 国产精品野战在线观看 | 嫁个100分男人电影在线观看| 丝袜美足系列| 久久精品aⅴ一区二区三区四区| 欧美人与性动交α欧美精品济南到| 美国免费a级毛片| 国产91精品成人一区二区三区| 久久久精品国产亚洲av高清涩受| 757午夜福利合集在线观看| 在线免费观看的www视频| 人人妻人人澡人人看| 国产精品国产高清国产av| 免费高清视频大片| 久久精品国产99精品国产亚洲性色 | 一级片免费观看大全| 久久国产亚洲av麻豆专区| 美女大奶头视频| 久久九九热精品免费| 老司机深夜福利视频在线观看| 一进一出好大好爽视频| 18禁黄网站禁片午夜丰满| 国产精品乱码一区二三区的特点 | 亚洲精品一二三| 精品久久久久久电影网| 国产成人精品久久二区二区免费| 午夜两性在线视频| 久久国产精品男人的天堂亚洲| 日本 av在线| 国产精品免费一区二区三区在线| 欧美午夜高清在线| 亚洲成人免费电影在线观看| 看黄色毛片网站| 激情在线观看视频在线高清| 在线观看一区二区三区| 一区二区三区精品91| 亚洲专区字幕在线| 亚洲熟妇中文字幕五十中出 | 国产av在哪里看| 午夜福利在线观看吧| 丝袜美足系列| 一a级毛片在线观看| 在线观看免费午夜福利视频| 又黄又爽又免费观看的视频| 日本黄色日本黄色录像| 人妻久久中文字幕网| 999精品在线视频| 亚洲七黄色美女视频| 国产99久久九九免费精品| 免费少妇av软件| 亚洲色图 男人天堂 中文字幕| 国产成+人综合+亚洲专区| 亚洲第一青青草原| 久久中文看片网| 一级作爱视频免费观看| 99国产精品一区二区三区| 亚洲,欧美精品.| 最新在线观看一区二区三区| 久久久久国产精品人妻aⅴ院| 不卡一级毛片| 久久中文字幕人妻熟女| 黄片播放在线免费| 亚洲欧美日韩无卡精品| 午夜福利免费观看在线| 免费一级毛片在线播放高清视频 | 亚洲国产欧美网| 国产亚洲精品第一综合不卡| 欧美精品一区二区免费开放| 757午夜福利合集在线观看| 十八禁网站免费在线| 久久久久久亚洲精品国产蜜桃av| 精品国产亚洲在线| 性色av乱码一区二区三区2| 91精品国产国语对白视频| 国产av一区二区精品久久| 亚洲成a人片在线一区二区| 一区二区三区激情视频| 91av网站免费观看| 精品国内亚洲2022精品成人| 日韩欧美一区视频在线观看| 日本wwww免费看| 国产精品亚洲一级av第二区| 国产精品av久久久久免费| 美女福利国产在线| 国产蜜桃级精品一区二区三区| 最新在线观看一区二区三区| 欧美一区二区精品小视频在线| av超薄肉色丝袜交足视频| 色综合婷婷激情| 18禁国产床啪视频网站| aaaaa片日本免费| 精品福利观看| 香蕉丝袜av| 热re99久久国产66热| 在线免费观看的www视频| 精品一区二区三卡| 在线观看午夜福利视频| 欧美乱码精品一区二区三区| 男女床上黄色一级片免费看| 国内久久婷婷六月综合欲色啪| 天天躁狠狠躁夜夜躁狠狠躁| av福利片在线| 亚洲精品在线美女| 久久人妻福利社区极品人妻图片| 久久人人爽av亚洲精品天堂| av欧美777| 成年女人毛片免费观看观看9| 91字幕亚洲| 一区在线观看完整版| 一级a爱视频在线免费观看| 香蕉丝袜av| 日日夜夜操网爽| av视频免费观看在线观看| 首页视频小说图片口味搜索| 中文亚洲av片在线观看爽| 日本a在线网址| 脱女人内裤的视频| 久久中文看片网| 伦理电影免费视频| 一区福利在线观看| 首页视频小说图片口味搜索| www.自偷自拍.com| 69精品国产乱码久久久| 欧美成人午夜精品| 免费人成视频x8x8入口观看| 无限看片的www在线观看| 国产av精品麻豆| 亚洲中文日韩欧美视频| 欧美激情极品国产一区二区三区| 最新美女视频免费是黄的| 亚洲国产精品999在线| aaaaa片日本免费| 自线自在国产av| 欧美日韩视频精品一区| x7x7x7水蜜桃| 精品国产超薄肉色丝袜足j| 高清毛片免费观看视频网站 | 亚洲欧美一区二区三区久久| av有码第一页| 黑人欧美特级aaaaaa片| 国产精品免费一区二区三区在线| 老司机午夜福利在线观看视频| 18禁观看日本| 午夜精品在线福利| 中文亚洲av片在线观看爽| cao死你这个sao货| 日韩国内少妇激情av| 国内毛片毛片毛片毛片毛片| 中文字幕高清在线视频| www.999成人在线观看| 国产aⅴ精品一区二区三区波| 婷婷六月久久综合丁香| 精品日产1卡2卡| 涩涩av久久男人的天堂| 免费在线观看黄色视频的| 黄网站色视频无遮挡免费观看| 免费看十八禁软件| 欧美精品亚洲一区二区| 在线看a的网站| 热99国产精品久久久久久7| 精品免费久久久久久久清纯| 日韩免费高清中文字幕av| av在线天堂中文字幕 | 侵犯人妻中文字幕一二三四区| 国产精品 国内视频| 水蜜桃什么品种好| 亚洲黑人精品在线| 激情在线观看视频在线高清| 久久中文字幕人妻熟女| 9191精品国产免费久久| 免费高清在线观看日韩| 91在线观看av| 久久久精品欧美日韩精品| 黄色毛片三级朝国网站| 国产欧美日韩精品亚洲av| 我的亚洲天堂| 亚洲精品美女久久久久99蜜臀| 精品一区二区三区av网在线观看| 国产精品久久久人人做人人爽| 色尼玛亚洲综合影院| 亚洲激情在线av| 中文字幕最新亚洲高清| 少妇被粗大的猛进出69影院| 日韩欧美在线二视频| 一区二区日韩欧美中文字幕| 多毛熟女@视频| 久久精品国产清高在天天线| 老司机午夜福利在线观看视频| 久久国产乱子伦精品免费另类| 亚洲免费av在线视频| 国产野战对白在线观看| 最新在线观看一区二区三区| 亚洲精品一区av在线观看| 性少妇av在线| 涩涩av久久男人的天堂| 久久久国产成人精品二区 | 无人区码免费观看不卡| 999久久久国产精品视频| 国产色视频综合| 黄色 视频免费看| 操出白浆在线播放| 大型黄色视频在线免费观看| 久久午夜亚洲精品久久| 欧美日韩亚洲高清精品| 人人妻人人澡人人看| 97碰自拍视频| 成人国语在线视频| 国产又爽黄色视频| 免费久久久久久久精品成人欧美视频| 久热这里只有精品99| 国产精品永久免费网站| 在线永久观看黄色视频| 制服人妻中文乱码| 免费在线观看完整版高清| 午夜免费鲁丝| 午夜福利免费观看在线| 咕卡用的链子| 久久人妻熟女aⅴ| 欧美不卡视频在线免费观看 | 欧美日韩福利视频一区二区| ponron亚洲| 精品熟女少妇八av免费久了| 国产男靠女视频免费网站| 自拍欧美九色日韩亚洲蝌蚪91| 精品一区二区三卡| 黄色丝袜av网址大全| 欧美黄色片欧美黄色片| av国产精品久久久久影院| 国产免费现黄频在线看| 一级毛片高清免费大全| 叶爱在线成人免费视频播放| 欧美国产精品va在线观看不卡| 久久久久久久精品吃奶| 亚洲欧美激情在线| 欧美人与性动交α欧美精品济南到| www.自偷自拍.com| 免费不卡黄色视频| 久久久国产欧美日韩av| 法律面前人人平等表现在哪些方面| 久久天堂一区二区三区四区| 人人澡人人妻人|