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

    超高應(yīng)變率下巖石的破壞機(jī)理

    2021-03-11 08:49趙高峰孫建華李世金張奔
    土木建筑與環(huán)境工程 2021年1期

    趙高峰 孫建華 李世金 張奔

    摘 要:通過對大理石的輕氣炮試驗進(jìn)行三維數(shù)值仿真重現(xiàn),探索了巖石在超高應(yīng)變率下的破壞機(jī)理。大理石試樣中兩個高速應(yīng)力計的實測應(yīng)力信號被用作數(shù)值模擬的匹配目標(biāo),結(jié)果表明,在高速沖擊壓縮載荷下,細(xì)觀靜水壓破壞是大理石的主要破壞機(jī)制,而不是通常認(rèn)為的細(xì)觀剪切或拉伸破壞。此外,為再現(xiàn)不同沖擊速度下大理石的物理實驗觀測數(shù)據(jù),提出了新的考慮時間非局部效應(yīng)的細(xì)觀速率相關(guān)巖石破壞模型,并且發(fā)現(xiàn)對巖石輕氣炮的數(shù)值模擬需采用準(zhǔn)靜態(tài)測量的彈性力學(xué)參數(shù)而不是超聲測量的彈性參數(shù)。最后對具有不同孔隙度的虛擬巖石樣本進(jìn)行進(jìn)一步的輕氣炮數(shù)值試驗,結(jié)果表明,巖石的孔隙度越大,在沖擊壓縮載荷下的動態(tài)強(qiáng)度和動力學(xué)率效應(yīng)越小。

    關(guān)鍵詞:動態(tài)破壞;輕氣槍試驗;沖擊壓縮荷載;格元模型

    1 Introduction

    The light-gas gun, first developed for hypervelocity research in aeronautics, was later used in the field of material science. The light-gas gun test is mainly used in material science for: 1) obtaining the Hugoniot curve of the material[1-3], 2) measuring dynamic failure strength[4], 3) studying the high-pressure phase transition[5], and 4) investigating impact-induced chemical reactions[6]. The light-gas gun test has also been applied to rock and rock-like materials. Grote et al.[7] conducted a flat plate impact test on cement mortar and concrete. They reported that the average flow stress of cement mortar and concrete increased to 1.3 GPa and 1.7 GPa, respectively, compared with the unconfined quasi-static compressive strength of 46 MPa and 30 MPa. The dynamic effect of reinforced concrete under shock compressive loading was also explored using the light-gas gun test[8]. The bearing capacity of the reinforced concrete was found to be improved with the increase of impact velocity and reinforcement ratio[8]. Recently, similar phenomena were observed by Zhang et al.[1] with marble and gabbro.

    The dynamic elasticity theory for solid materials was adopted to interpret experimental data of the light-gas gun test[1,7-8]. The constitutive model based on the dynamic elasticity theory can be used with the finite element method (FEM) for numerical simulation of the light-gas gun test. Liu et al.[9] simulated the dynamic fracture mode of tungsten alloy by LS-DYNA (an explicit FEM code) and captured the characteristics in their experimental results. Lopatnikov et al.[10] simulated the light-gas gun test of foamed aluminum plates using LS-DYNA and demonstrated that the dynamic deformation-time relationship and kinetic energy changes of the impact plates at different impact velocities predicted by the FEM were consistent with their theoretical model. Nevertheless, as pointed out by Riedel et al.[11], the intrinsic heterogeneity of rock-like material might result in challenges to the determination of the equation of state parameters and the need for a mesomechanical model. Duan et al.[12] incorporated a statistical isotropic elastic micro-crack model into the FEM to simulate the propagation of plane shock waves in soda-lime glass and reproduced their experimental data. Recently, besides FEM with the mesomechanical model[13-14], a number of discontinuum-based numerical methods[15-20] have also been applied in the study of the dynamic failure of rock. The molecular dynamics (MD)[15] and the lattice type model[18-20] were successfully applied to study the dynamic failure of rock. However, most of these studies only reproduced the failure mode of the rock using a mesoscopic constitutive model with tensile failure. The implementation and practicability of the model to simulate the light-gas gun test is still unclear.

    In this work, we used a lattice type model to explore the light-gas gun test conducted on marble by Zhang et al.[1]. Our main purposes were to answer the following questions:

    1) What kind of mesoscopic damage information can be obtained from the experimental data of the light-gas gun test?

    2) What kind of elastic parameters, the dynamic elastic ones or the quasi-static elastic ones, control the shock loading propagation within the rock specimen?

    3) How to describe the rate dependency of rock under shock compressive loading in the light-gas gun test using a mesoscopic constitutive model?

    4) What is the role of the rock's mesostructure in the dynamic strength and dynamic effects (strain/loading rate dependency) of rock under shock compressive loading?

    In this work, the distinct lattice spring model (DLSM) was adopted as the numerical tool due to its advantages in modeling rock failure based on simple mechanical elements, e.g., Newton's second law and springs with simple constitutive models. These characteristics of the DLSM make it a suitable numerical tool for failure mechanism study. In addition, the DLSM is computationally appropriate for full 3D simulation. The following sections of this paper include a brief introduction of the physical test and numerical model, a brief description of the development of rate-dependent constitutive models and a comparison of model response with published data of the light-gas gun test conducted by Zhang et al.[1], based on which the failure mechanism of rock under shock compressive loading is explored. Finally, a few conclusions are derived to answer the aforementioned questions.

    2 Methods

    2.1 The light-gas test

    The light-gas test data reported by Zhang et al.[1], obtained using the experimental facilities at the Cavendish Laboratory, Cambridge, UK[21], are adopted in this work. Fig.1 shows the basic working principle of the experimental facilities as well as the components of the marble specimen and the copper flyer. Different from the traditional rock mechanics tests, two stress gauges were placed in the composite specimen made of marble and PMMA (see Fig.1(b)). The two stress gauges were used to obtain the stress history curve of the corresponding location inside the specimen during the test, that is, the waveform of the shock wave. Since there is a certain distance between the two stress gauges, the velocity of the shock wave in the rock specimen can be obtained based on these signals. After the light-gas gun test, both the composite rock specimen and the flyer became powdery and dissipated, as a result, no morphological data was recorded for the failure process, or the failure pattern of the rock specimen. Another important parameter of the light-gas gun test is the impact velocity of the copper flyer, which is the only controllable input of the light-gas gun test. A number of different impact velocities were adopted by Zhang et al.[1] to obtain the Hugoniot parameters of the marble. Basic material properties of the marble, PMMA and copper flyer were also provided in reference [1]. The longitudinal and shear wave velocities were obtained using an ultrasonic transducer, which can be further used to obtain the dynamic elastic parameters. The quasi-static elastic parameters of the marble were obtained by the standard uniaxial compression test, the parameters of which are listed in Table 1. The dynamic elastic parameters of PMMA and copper were close to their quasi-static counterparts, therefore, no quasi-static tests were conducted on the copper and PMMA.

    2.2 Distinct Lattice Spring Model (DLSM)

    The lattice spring model(LSM) was first developed by Hrennikoff in 1941[22]. Its basic principle is to represent the mechanical responses of a solid through a group of spring-like interactions. The LSM is regarded as the ancestor of both the FEM and the discrete element method (DEM). Due to its simplicity, the LSM has been widely used in the study of many fundamental mechanical phenomena of solids[23-25]. The distinct lattice spring model (DLSM) was developed by Zhao et al.[26] to overcome the Poisson's limitation in the classical LSM. As shown in Fig.2(a), the basic principle of the DLSM is to represent the solid as a group of particles linked through spring bonds, which consist of a normal spring and a shear spring. Defining the normal unit vector: n=(nx ny nz)T is directed from particle i to particle j, which are connected by a normal spring, and the normal deformation of the spring is defined as

    in which ks is the stiffness of the shear spring and u*s is the ultimate shear deformation.

    Eqs.(2) and (4) are the constitutive models adopted in the DLSM for brittle solids. Fig. 3(a) and (b) show the constitutive model for the normal spring and the constitutive model for the shear spring, respectively. Since a solid is represented as a spring network in the DLSM (as shown in Fig. 2(a)), the fracturing process is manifested by the gradual breakage of these springs (mesoscopic failure events). It has the following advantages in simulating the fracturing of brittle materials. First, these constitutive models are easily implemented, and the simulation results are easily explained, as fewer parameters are required compared to other, discontinuous models, e.g., the bonded DEM. Third, the straight forward parallelization of the DLSM[27-28] gives the model relatively high computational performance.

    Another distinctive feature of the DLSM is that the relationship between the spring parameters(kn and ks) and the macroscopic material properties (E and v) is derived from the hyper-elastic theory. No calibration is required for the mesoscopic elastic parameters, as they can be calculated using the following equations.

    where li is the initial length of the ith spring and Vm is the volume representing the computational model. More details of the DLSM and its recent development can be found in reference [25-27].

    2.3 Computational model for the light-gas gun test

    Fig.2(b) shows the computational model for the light-gas gun test conducted by Zhang et al.[1] on a marble specimen. Here, the marble specimen is represented as a group of mass points linked by springs. The dynamic contact between the particles was considered in the DLSM, so the impacts on the light-gas gun test could be directly simulated. The influences of thermal and chemical reactions were not considered in this work. Therefore, its mechanical response is characterized by a mass-spring system based on Newton′s second law and Hooke's law[26]. As shown in Fig.2(b), the computational model consists of a copper flyer, a composite marble specimen and a PMMA plate. This three-dimensional numerical model can be viewed as a numerical test parallel to Zhang et al[1], as shown in Fig.1(b). The specific geometric parameters are consistent with the physical experiment. The only boundary condition input is the corresponding impact velocity of the copper flyer. For a closer comparison with the outputs of the actual experiment in Zhang et al.[1], combined with the characteristics of the DLSM, a numerical stress gauge scheme was developed to record the stress data in the simulation. As shown in Fig.2(c), a rectangular plane was placed in the model, and its normal vector is nA. During the calculation, assume the normal unit direction vector of the spring bond cut by this plane is nbi, and the total spring bond force accounting for the normal and shear springs is Fbi(t)=Fni(t)+Fsi(t), then the corresponding stress of the numerical stress gauge is given as

    Eq. (8) gives the time history of the normal stress on the stress gauge.

    The particle size used in this work is 1 mm and there are about 100 000 particles. The input parameters of the light-gas gun test are based on experimental tests. The output of the numerical test is controlled by the selection of the meso-mechanical constitutive model and the corresponding constitutive parameters. A key point of this work is to minimize the difference between the numerical prediction and the physical experimental result. Thus, the failure mechanism of the rock under shock compressive loading may be interpreted from the searching process of the meso-mechanical constitutive model and its parameters.

    2.4 Tri-linear hydro-compressive model

    Many researchers suggest that macroscopic compression and shear failure is essentially tensile failure at the mesoscopic level[28-30]. However, most of those studies are limited to tensile failure of the normal spring[18, 20, 30]. Fig.3(a) shows the simple brittle constitutive model of the normal spring considering only tensile failure. More comprehensive forms, e.g., considering the damage and plastic deformation, can be developed. Jiang and Zhao[31] recently developed a coupled damage plasticity model to describe the dynamic crack propagation of a gypsum-like 3D printing material. When considering the interaction between two particles, shear failure is a natural logical extension. However, there are very few studies on shear failure in classical LSMs due to the absence of shear interaction. In the DLSM, shear failure can be considered due to the introduction of the multi-body shear spring. The corresponding constitutive model is illustrated in Fig.3(b). Shear failure was widely adopted in the bonded DEM, which usually takes into account both the influence of the cohesion of the contact bond between two particles and the frictional angle. The shear constitutive model presented in Fig.3(b) can be viewed as a special case when the frictional angle is zero. Therefore, the prediction of shear fracturing of the DLSM using the constitutive model shown in Fig.3(b) is not conservative.

    When considering the interaction between two particles, it can be suggested that there is another possibility of damage, that is, the compression failure of the normal spring, which is related to the failure of the material under hydrostatic compression. If there is no hydrostatic failure involved, there is no need to introduce the hydro-compressive failure model into the DLSM. In this work, a trilinear constitutive model (see Fig.4) is introduced to describe the hydrostatic failure of the rock. Four parameters are required to determine the constitutive model without considering the loading and unloading. The first parameter, uc1, represents the first hydro-compression yield of the normal spring. The hardening or softening stage starts, followed by the normal spring yields. The second parameter, βredc1, corresponds to the reduction factor of the spring stiffness of the first yield. When it is greater than zero, it represents hardening, and when it is smaller than zero, it represents softening. The other two parameters, uc2 and βredc2, correspond to the secondary hydro-compression yield deformation and the stiffness reduction coefficient, respectively. The mathematic equation of this tri-linear constitutive model is as follows.

    The response of the constitutive model is controlled by these four parameters,uc1,uc2,βredc1and βredc2

    To introduce the effects of plasticity and damage, the same approach as that adopted by Jiang and Zhao[31] is used. First, rewrite Eq.(10) as follows.

    where S(u) is a state parameter of the normal spring, which is defined as follows.

    The interpretation of Eq.(10) can be expressed as the following formula if it is based on damage mechanics.

    Then, the damage function can be expressed as follows.

    In the calculation of the spring force, the maximum damage variable corresponding to the loading history is recorded as D*(u), then the damage constitutive equation considering the loading and unloading can be written as

    The loading and unloading response of this constitutive model is as shown in Fig.4(a). If Eq.(10) is regarded as a plastic response, it can be expressed as

    Compare Eq.(14), we will have the plastic deformation up as

    Similarly, to record the maximum plastic deformation up* experienced by the spring at the moment of loading, the pure plastic constitutive model considering loading and unloading can be written as

    The corresponding pure plastic constitutive model is shown in Fig.4(b). Combined with Eq.(15) and (18), a coupled damage-plastic model can be constructed by introducing the damage-plastic coupling coefficient as follows.

    When λdp=1, it represents the pure damage constitutive model, and when λdp=0, it represents the pure plastic constitutive model. In this work, we consider the copper behavior as pure plastic material, while the rock has a pure damage response.

    Under dynamic loading, the material has obvious dynamic rate effects. To capture this phenomenon, it is convenient to introduce a rate-dependent model. The approach developed by Zhao et al.[18] is adopted in this work to bring rate dependency into the tri-linear constitutive model. The idea is to change the compressive strength parameters of the trilinear constitutive model as a function of the spring deformation rate. It was found that the instantaneous value of the spring deformation rate may not be able to reproduce the correct macroscopic dynamic effect[18]. To solve this problem, Zhao et al.[18] proposed the concept of time non-localization, which uses the average value of the spring deformation rate from the start time to the current time. In this work, following a similar idea, the local average deformation rate is given as

    The average deformation rate of the spring is only calculated when the spring is deformed beyond 0.99uc1. Based on this concept, the strength parameters of the corresponding dynamic tri-linear constitutive model are given by the following formula.

    where udync1 is the first yield deformation considering the dynamic effect, udync2 is the second yield deformation considering the dynamic effect, vref is the reference deformation rate, and α is a dimensionless coefficient. Fig.5(a) shows three different loading curves; each represents a different loading rate (deformation rate). The loading rate has alternating positive and negative values, which is common in dynamic loads. Fig.5(b) shows the dynamic response based on the proposed dynamic tri-linear constitutive model, which shows good robustness for dynamic/cyclic loading. The dynamic constitutive model contains a total of six parameters with clear geometric and physical meanings. The selection of the parameters in the actual simulation process will be further illustrated in the following section.

    3 Numerical modeling and discussion

    3.1 Elastic parameters selection

    A proper selection of parameters is fundamental to numerical simulation. For the light-gas gun test, Zhang et al.[1] performed two kinds of test to obtain the elastic parameters of their marble specimens. They conducted ultrasonic tests and obtained the wave velocities, which is a main variable of the dynamic elastic modulus and Poisson's ratio, as shown in Table 1. The quasi-static elastic modulus and Poisson's ratio of the marble specimens were also obtained through the traditional uniaxial compression test (see Table 1). The quasi-static elastic parameters of the copper and PMMA are close to the dynamic ones, therefore, they were not tested and reported in reference [1]. Now, the challenge is the selection of the elastic parameters, dynamic or quasi-static ones, for the numerical modeling of the light-gas gun test. Here, the light-gas gun test with an impact velocity of 385 m/s is used as our modeling target. During the test, stress signals were recorded by two stress gauges, which were used to obtain the characteristics of the shock wave propagation within the specimen. It is reasonable to use the elastic parameters for ultrasonic wave propagation within the specimen to simulate the shock wave propagation. Therefore, dynamic elastic parameters are assigned to the marble as the first step, i.e., E=55 GPa, v=0.31, for the numerical simulation. Here, a pure elastic dynamic simulation is performed by setting all failure parameters extremely large to prevent any damage from occurring, e.g. u*n=1 mm and u*s=1 mm. In order to achieve a quantitative comparison between the experimental data and the numerical data, an adjustment of the initial recording time is performed to make the initial waveform of the first numerical stress gauge match the corresponding physical signal. The time shift value was also applied to the signal of the second numerical stress gauge. As shown in Fig.6(a), the initial stage of the numerical prediction of the second numerical stress gauge is significantly earlier than the experimental observation. This indicates that the numerically predicted shock wave velocity is markedly larger than that of the physical observation when the dynamic elastic parameters are used for the light-gas gun test. The quasi-static elastic parameters for the marble specimen, i.e., E=40 GPa, v=0.20, are also adopted. The simulation results are shown in Fig.7(b). It is noted that a better fitting is observed. Therefore, in the subsequent numerical modeling, for the marble specimen, the quasi-static elastic parameters are selected as the input values. In terms of the stress peaks, the numerical results are also much larger than the experimental observation, revealing that damage/failure occurred within the specimen during the light-gas gun test, otherwise a much higher peak value would be detected in the physical test. In the following section, we will further explore the corresponding failure mechanism of the specimen under shock compressive loading.

    3.2 Failure mechanism of the marble specimen in the light-gas gun test

    In this section, the mesoscopic failure mechanism of the marble specimen during the light-gas gun test conducted by Zhang et al.[1] is explored. First, the failure parameters were set as u*n=0.000 1 mm, and u*s=1 mm, which means only tensile failure was considered. According to the empirical equation proposed by Zhao[32], the represented macro tensile strength is about u*nE/d-=4 MPa, where E is the elastic modulus and d- is the mean particle size of the computational model. Fig.7(a) shows the corresponding simulation results. Fig.8 shows the shock wave propagation and damage evolution of the marble specimen when considering only the mesoscopic tensile failure. Initially, the copper sheet collides with the rock sample at a high speed to generate a shock wave, and contact occurs at 3.00 microseconds. As shown in Fig.8, the speed of the copper sheet is immediately transmitted to the rock specimen and the part in contact with the copper sheet immediately breaks down. With the propagation of the shock wave, the damaged area spreads continuously around the copper piece (Fig.8(b)). The damage first propagates in the direction of impact, then it expands radially, and finally, the entire specimen is damaged. However, according to the observation presented in Fig.7(a), the influence on the stress waveforms at the last stage is minor. The amplitude of the stress wave recorded at the first stress gauge is not reduced at all, even if the marble specimen is fully damaged in tensile failure. Considering the mesoscopic tensile failure has little effect on the experimental observations, it is not attributable to the failure mechanism of the marble specimen in the light-gas gun test. It should be mentioned that using the DLSM the uniaxial compressive failure can be reproduced with only mesoscopic tensile failure[32]. It is speculated that the failure mechanism of the marble specimen may differ from the compressive failure observed in the classical uniaxial compression test. To further investigate the influence of shear failure between particles, the failure parameters are set as u*n=1 mm, and u*s=0.000 1 mm. The corresponding numerical calculation results are shown in Fig.7(b). In the DLSM, the shear constitutive model corresponds to the case where the friction coefficient between particles is zero, that is, the shear failure is fully considered. From the observations of Fig.7(b), it is concluded that the mesoscopic shear failure is also not able to account for the failure mechanism.

    In the light-gas gun test, the copper flyer was completely destroyed. Whether the recorded waveform can be affected by the damage to the copper flyer needs to be investigated. Here one may consider the compressive failure of the normal spring of the copper using the hydro-compressive constitutive model, as shown in Fig.4. The corresponding failure parameters were selected as uc1=0.001 mm, uc2=0.200 mm, βredc1=0, βredc2=0 and λdp=0 to represent the ideal elastoplastic response of the copper. Fig.9(a) shows the numerical simulation results. Compared with the experimental results, the peak is still unable to meet the experimental value. The compressive failure parameters of the marble then were taken as uc1=0.003 mm, uc2=0.200 mm, βredc1=0, βredc2=0 and λdp=1. The compression failure parameter uc1=0.003 mm, corresponds to the macroscopic failure strength of about 120 MPa, which is the typical value of the uniaxial compression strength of marble. Fig.9(b) shows the corresponding numerical simulation results. By comparing the numerical and physical results, it is concluded that the compression failure between particles in the light-gas gun test is caused by the hydro-static compression failure of the rock rather than uniaxial compression failure. In the following section, the compressive constitutive parameters will be adjusted to make the computational model reproduce the experimental results.

    The tri-linear constitutive model developed in this work consists of four parameters. First, the values of these four parameters are taken as uc2=0.200 mm, βredc1=0,βredc2=0 and λdp=1. By continuously adjusting uc1and performing numerical modeling, it is found that when uc1 is 0.030 mm, the peak value of the numerical shock wave of stress gauge 1 is closest to the experimental results (see Fig.10(a)). However, the reduction of the stress wave of the second stress gauge was obvious. This attenuation is caused by the relatively coarse particle size adopted in the simulation. For a full three-dimensional model, there are already 110 000 particles at a particle size of 1 mm. Considering the computing capacity limitation of the available hardware and that the main purpose is to explore the failure mechanism rather than precisely capture the attenuation of the shock wave, the computational model with the particle size of 1 mm is preferred in the following simulation with a focus on the comparison of the stress waveform at stress gauge 1. Besides the attenuation, the platform segment of the stress wave at stress gauge 1 was also not well reproduced. By further adjusting other parameters, it is found that by setting βredc1=0.1, uc2=0.1 mm, and βredc2=0, the numerical simulation results give a closer fitting to the physical test (see Fig.10(b)), which indicates that a light hardening has occurred for the hydro-compression failure of the marble under shock compressive loading.

    3.3 Dynamic effect of the mesoscopic hydro-compressive failure

    So far,in the experimental data of Zhang et al.[1] while the impact velocity of 385 m/s is well simulated, one needs to generalize it for other impact velocities adopted by Zhang et al. The strength Sx can be determined as the peak value of the stress waveform recorded at numerical stress gauge 1. Fig.11 shows the experimental data of the dynamic strength of the marble Sx versus different particle velocities. As suggested by Zhang et al.[1], the relationship between impact velocity vi and particle velocity vp may be estimated as vp=0.67vi, which is adopted in the data process of this simulation. Additional numerical simulations using different impact velocities, i.e., 550 m/s, 650 m/s, 800 m/s and 1 000 m/s, were conducted using the same computational model and constitutive parameters. The rate-independent (RI) constitutive model was first used. As shown in Fig. 11, the numerical prediction using the DLSM with the RI constitutive model failed to reproduce the dynamic effect of the strength value observed in the light-gas gun test. Then, the rate dependent (RD) constitutive model is used, in which the strength parameters of the tri-linear constitutive model are described by Eq.(21) and (22). Fig.11 shows the numerical results when v0ref=255 mm/s and α=0.43. From the numerical modeling results, it is concluded that the mesoscopic hydro-compressive constitutive model has to include the rate-dependent effect. In other words, the mesoscopic failure of the marble under shock compressive loading includes a dynamic effect.

    3.4 Influence of mesostructure

    In this section, the influence of the mesostructure on the dynamic failure of rock under shock loading is explored. Fig.12 shows the corresponding computational models for the rock specimens with different mesostructures. These computational models are formed by randomly removing a given number of particles from the marble specimen. The porosity, which was used to characterizs different computational models, is defined as n=number of particles being excavated / total number of particles in the original model. As shown in Fig.12, four pore models with a porosity of 5%, 10%, 15%, and 20% are constructed. All the constitutive parameters are the same as those of the previous section. Fig.13 shows the numerical results of the predicted strength of those different porosity models under different impact velocities. Overall, larger porosity results in lower dynamic strength, which seems logically correct. However, there are some variations, for example, when porosity is 5%, a higher strength is obtained. The reason might be that the inter-particle velocities become more violent under shock compressive loading when a small number of particles were removed from the original model. It might trigger the dynamic effect of the constitutive model and consequently, result in a higher strength. Nevertheless, when the porosity continuously increases, the inter-particle velocity increase may be released in the lateral direction due to the existence of a large mesoscopic free surface. This may result in the increase of strength due to the dynamic effect becoming inconspicuous (as shown in Fig.13).

    4 Conclusions

    In this work,light-gas gun tests on marble were numerically investigated using the DLSM to study the failure mechanism of rock under compressive shock loading. Based on a detailed comparison of the numerical response with published experimental data, it was found that the elastic parameters controlling the compressive shock wave propagation in the rock specimen should be determined through quasi-static tests rather than ultrasonic tests. Moreover, mesoscopic tensile failure and shear failure between rock grains are not the mechanisms of rock failure under compressive shock loading. To reproduce the shock stress waveform observed in the light-gas gun test, a mesoscopic hydro-compressive constitutive model is needed. In other words, it is reasonable to conclude that the light-gas gun test may be able to provide information about the hydro-compressive failure of rock under dynamic loading and provide calibration data for the DLSM or other numerical methods to determine the constitutive equations under dynamic hydro-compression loading conditions. It is noted that a rate-dependent hydro-compressive constitutive model is also necessary to reproduce the experimental observation of the strength increase under different impact velocities. Finally, the influence of the mesostructure in terms of porosity was found to decrease both the dynamic strength and the dynamic effect. The findings in this work may provide a better understanding of rock failure under compressive shock loading, which might be useful for a more rational protective design of underground rock engineering structures under blasting loading[33-36].

    Acknowledgements

    This research is financially supported by the National Key R&D Program of China (Grant No. 2018YFC0406804) and the National Natural Science Foundation of China (Grant No. 11772221).References:

    [1] ZHANG Q B, BRAITHWAITE C H, ZHAO J. Hugoniot equation of state of rock materials under shock compression [J]. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2017, 375(2085): 20160169.

    [2] ROSENBERG Z, YAZIV D, YESHURUN Y, et al. Shear strength of shock-loaded alumina as determined with longitudinal and transverse manganin gauges [J]. Journal of Applied Physics, 1987, 62(3): 1120-1122.

    [3] GEBBEKEN N, GREULICH S, PIETZSCH A. Hugoniot properties for concrete determined by full-scale detonation experiments and flyer-plate-impact tests [J]. International Journal of Impact Engineering, 2006, 32(12): 2017-2031.

    [4] GRADY D E. The spall strength of condensed matter [J]. Journal of the Mechanics and Physics of Solids, 1988, 36(3): 353-384.

    [5] DUNN J E, FOSDICK R, SLEMROD M. Shock induced transitions and phase structures in general media [M]. Springer, 1993.

    [6] SEKINE T. Shock wave chemical synthesis [J]. European Journal of Solid State and Inorganic Chemistry, 1997, 34: 823-833.

    [7] GROTE D L, PARK S W, ZHOU M. Dynamic behavior of concrete at high strain rates and pressures: I. experimental characterization [J]. International Journal of Impact Engineering, 2001, 25(9): 869-886.

    [8] LIU H F, FU J, NING J G. Experimental study on the dynamic mechanical properties of reinforced concrete under shock loading [J]. Acta Mechanica Solida Sinica, 2016, 29(1): 22-30.

    [9] LIU H Y, SONG W D, REN H L. Dynamic response of tungsten-nickel-iron composites under impact loadings [J]. International Journal of Nonlinear Sciences and Numerical Simulation, 2009, 10(8): 993-1004.

    [10] LOPATNIKOV S L, GAMA B A, HAQUE M J, et al. High-velocity plate impact of metal foams [J]. International Journal of Impact Engineering, 2004, 30(4): 421-445.

    [11] RIEDEL W, WICKLEIN M, THOMA K. Shock properties of conventional and high strength concrete: Experimental and mesomechanical analysis [J]. International Journal of Impact Engineering, 2008, 35(3): 155-171.

    [12] DUAN Z P, ZHANG Y G, ZHANG L S, et al. Multi-thickness target plate impact experimental approach to failure waves in soda-lime glass and its numerical simulation [J]. International Journal of Nonlinear Sciences and Numerical Simulation, 2012, 13(1): 3-11.

    [13] TANG C A. Numerical simulation of progressive rock failure and associated seismicity [J]. International Journal of Rock Mechanics and Mining Sciences, 1997, 34(2): 249-261.

    [14] ZHU W C, LIU J, YANG T H, et al. Effects of local rock heterogeneities on the hydromechanics of fractured rocks using a digital-image-based technique [J]. International Journal of Rock Mechanics and Mining Sciences, 2006, 43(8): 1182-1199.

    [15] KRIVTSOV A M. Relation between spall strength and mesoparticle velocity dispersion [J]. International Journal of Impact Engineering, 1999, 23(1): 477-487.

    [16] SILLING S A. Reformulation of elasticity theory for discontinuities and long-range forces [J]. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175-209.

    [17] ZHOU X P, WANG Y T, QIAN Q H. Numerical simulation of crack curving and branching in brittle materials under dynamic loads using the extended non-ordinary state-based peridynamics [J]. European Journal of Mechanics-A, 2016, 60: 277-299.

    [18] ZHAO G F, RUSSELL A R, ZHAO X B, et al. Strain rate dependency of uniaxial tensile strength in Gosford sandstone by the Distinct Lattice Spring Model with X-ray micro CT [J]. International Journal of Solids and Structures, 2014, 51(7/8): 1587-1600.

    [19] ZHANG Z N, YAO Y, MAO X B. Modeling wave propagation induced fracture in rock with correlated lattice bond cell [J]. International Journal of Rock Mechanics and Mining Sciences, 2015, 78: 262-270.

    [20] ZHAO G F, XIA K W. A study of mode-I self-similar dynamic crack propagation using a lattice spring model [J]. Computers and Geotechnics, 2018, 96: 215-225.

    [21] BOURNE N K, ROSENBERG Z, JOHNSON D J, et al. Design and construction of the UK plate impact facility [J]. Measurement Science and Technology, 1995, 6(10): 1462-1470.

    [22] HRENNIKOFF A. Solution of problems of elasticity by the framework method [J]. Journal of Applied Mechanics, 1941, 8: A619-A715.

    [23] FINEBERG J, GROSS S P, MARDER M, et al. Instability in dynamic fracture [J]. Physical Review Letters, 1991, 67(4): 457.

    [24] GUILLARD F, GOLSHAN P, SHEN L M, et al. Dynamic patterns of compaction in brittle porous media [J]. Nature Physics, 2015, 11(10): 835-838.

    [25] ZHAO G F. Developing a four-dimensional lattice spring model for mechanical responses of solids [J]. Computer Methods in Applied Mechanics and Engineering, 2017, 315: 881-895.

    [26] ZHAO G F, FANG J N, ZHAO J. A 3D distinct lattice spring model for elasticity and dynamic failure [J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2011, 35(8): 859-885.

    [27] ZHAO G F, FANG J N, SUN L, et al. Parallelization of the distinct lattice spring model [J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2013, 37(1): 51-74.

    [28] ZHAO G F, CHEN H. Performance of the parallel four-dimensional lattice spring model using Alibaba cloud [J]. Journal of Civil and Environmental Engineering, 2019, 41(3): 1-11.

    [29] POTYONDY D O, CUNDALL P A. A bonded-particle model for rock [J]. International Journal of Rock Mechanics and Mining Sciences, 2004, 41(8): 1329-1364.

    [30] JIANG C, ZHAO G F. On crack propagation in brittle material using the distinct lattice spring model [J]. International Journal of Solids and Structures, 2017, 118/119: 41-57.

    [31] JIANG C, ZHAO G F. Implementation of a coupled plastic damage distinct lattice spring model for dynamic crack propagation in geomaterials [J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2018, 42(4): 674-693.

    [32] ZHAO G F. Development of micro-macro continuum-discontinuum coupled numerical method [D]. EPFL, Switzerland, 2010.

    [33] XIE H, ZHU J, ZHOU T, et al. Conceptualization and preliminary study of engineering disturbed rock dynamics [J]. Geomechanics and Geophysics for Geo-Energy and Geo-Resources, 2020. https://doi.org/10.1007/s40948-020-00157-x.

    [34] ZHOU X P, YANG H Q. Micromechanical modeling of dynamic compressive responses of mesoscopic heterogenous brittle rock [J]. Theoretical and Applied Fracture Mechanics, 2007, 48(1): 1-20.

    [35] LIU P F, ZHOU X P, QIAN Q H, et al. Dynamic splitting tensile properties of concrete and cement mortar [J]. Fatigue & Fracture of Engineering Materials & Structures, 2020, 43(4): 757-770.

    [36] ZHANG B, ZHAO G F. Dynamic failure and energy dissipation of rock ring specimen based on 4D lattice spring model [J]. Journal of Civil and Environmental Engineering, 2019, 41(2): 20-28.

    (編輯 章潤紅)

    av国产免费在线观看| 一个人观看的视频www高清免费观看| 亚洲欧美中文字幕日韩二区| 久久久久久久亚洲中文字幕| 插逼视频在线观看| 久久精品国产自在天天线| 久久久久久久午夜电影| 在线免费观看不下载黄p国产| 亚洲av成人精品一区久久| 夫妻性生交免费视频一级片| 成人国产麻豆网| 在线观看美女被高潮喷水网站| 国产精品熟女久久久久浪| 高清在线视频一区二区三区 | 一级黄色大片毛片| 亚洲av.av天堂| 岛国在线免费视频观看| 日本三级黄在线观看| 国产精品乱码一区二三区的特点| 久久精品久久久久久噜噜老黄 | 亚洲成人av在线免费| 成人三级黄色视频| 一个人观看的视频www高清免费观看| 如何舔出高潮| 国产黄片美女视频| 97人妻精品一区二区三区麻豆| 乱人视频在线观看| 国产一级毛片在线| 黄色一级大片看看| 26uuu在线亚洲综合色| 少妇裸体淫交视频免费看高清| 韩国av在线不卡| 天堂av国产一区二区熟女人妻| 欧美精品一区二区大全| a级毛片免费高清观看在线播放| 国产老妇伦熟女老妇高清| 国产 一区 欧美 日韩| 亚洲成人av在线免费| 精品一区二区免费观看| 国产精品久久久久久av不卡| 国产 一区 欧美 日韩| 九九在线视频观看精品| 国产成人91sexporn| 毛片一级片免费看久久久久| 亚洲av不卡在线观看| 精品少妇黑人巨大在线播放 | 又爽又黄a免费视频| 亚洲成色77777| 亚洲欧洲国产日韩| 最近手机中文字幕大全| 男女那种视频在线观看| 亚洲精品aⅴ在线观看| 国产成年人精品一区二区| 亚洲精品国产av成人精品| 级片在线观看| 精品熟女少妇av免费看| 国产精品三级大全| 成年免费大片在线观看| 高清视频免费观看一区二区 | 日本爱情动作片www.在线观看| 色视频www国产| 亚洲无线观看免费| 久久精品国产自在天天线| 91久久精品国产一区二区成人| 一卡2卡三卡四卡精品乱码亚洲| 国产探花极品一区二区| eeuss影院久久| 最近中文字幕高清免费大全6| 国产精品伦人一区二区| 九草在线视频观看| 亚洲自偷自拍三级| 热99re8久久精品国产| 欧美变态另类bdsm刘玥| 亚洲欧美日韩无卡精品| av卡一久久| av线在线观看网站| 亚洲成人av在线免费| 丰满人妻一区二区三区视频av| 一本久久精品| 99久久精品一区二区三区| 国产精品综合久久久久久久免费| 久久精品综合一区二区三区| 婷婷色av中文字幕| 亚洲av二区三区四区| 亚洲国产精品专区欧美| 插阴视频在线观看视频| 国产激情偷乱视频一区二区| 欧美日韩在线观看h| 一个人免费在线观看电影| 国产精品伦人一区二区| 国产老妇女一区| 菩萨蛮人人尽说江南好唐韦庄 | 岛国毛片在线播放| 国产私拍福利视频在线观看| 欧美色视频一区免费| 啦啦啦韩国在线观看视频| 国产乱人视频| 岛国在线免费视频观看| 国产成人精品一,二区| 大话2 男鬼变身卡| 欧美一区二区国产精品久久精品| 三级国产精品片| 亚洲,欧美,日韩| 亚洲精品成人久久久久久| 99国产精品一区二区蜜桃av| 国产精品麻豆人妻色哟哟久久 | 床上黄色一级片| 三级毛片av免费| 91精品一卡2卡3卡4卡| 国产真实伦视频高清在线观看| 最近中文字幕高清免费大全6| 蜜桃久久精品国产亚洲av| 老司机影院毛片| av又黄又爽大尺度在线免费看 | 国产av码专区亚洲av| 爱豆传媒免费全集在线观看| 麻豆乱淫一区二区| 国产精品1区2区在线观看.| 日韩一本色道免费dvd| 久久久久久大精品| 国产91av在线免费观看| 国产老妇女一区| 国产成人精品婷婷| 人妻系列 视频| 欧美日韩国产亚洲二区| 久久精品夜夜夜夜夜久久蜜豆| 97热精品久久久久久| 99热这里只有是精品50| 91aial.com中文字幕在线观看| 一区二区三区四区激情视频| 最近最新中文字幕大全电影3| 三级国产精品欧美在线观看| 午夜亚洲福利在线播放| 亚洲不卡免费看| 你懂的网址亚洲精品在线观看 | 成人高潮视频无遮挡免费网站| 热99在线观看视频| 国内精品宾馆在线| 亚洲,欧美,日韩| 亚洲怡红院男人天堂| 毛片女人毛片| 97超碰精品成人国产| 免费观看a级毛片全部| 亚洲五月天丁香| 日本一本二区三区精品| 欧美日本亚洲视频在线播放| 欧美成人一区二区免费高清观看| 欧美高清成人免费视频www| 国产精品无大码| 少妇猛男粗大的猛烈进出视频 | 高清毛片免费看| av.在线天堂| 亚洲国产欧美人成| 欧美高清性xxxxhd video| av.在线天堂| 日韩强制内射视频| 亚洲精品乱码久久久v下载方式| 高清午夜精品一区二区三区| 狂野欧美激情性xxxx在线观看| 蜜桃亚洲精品一区二区三区| 身体一侧抽搐| av视频在线观看入口| 亚洲天堂国产精品一区在线| 久久热精品热| 天天躁日日操中文字幕| 黄色欧美视频在线观看| 欧美bdsm另类| 日本欧美国产在线视频| 日本熟妇午夜| 日本猛色少妇xxxxx猛交久久| 亚洲一区高清亚洲精品| 亚洲中文字幕一区二区三区有码在线看| 久久久精品94久久精品| 久久精品国产鲁丝片午夜精品| 边亲边吃奶的免费视频| 久久99热这里只频精品6学生 | 在线观看一区二区三区| 狂野欧美白嫩少妇大欣赏| 日韩,欧美,国产一区二区三区 | 亚洲精华国产精华液的使用体验| 免费av毛片视频| 国产美女午夜福利| 国语自产精品视频在线第100页| 国产亚洲91精品色在线| 99久久人妻综合| 欧美精品一区二区大全| 成人午夜高清在线视频| 精品久久久噜噜| 村上凉子中文字幕在线| 亚洲欧美中文字幕日韩二区| 亚洲色图av天堂| 国产亚洲91精品色在线| 国产精品永久免费网站| 自拍偷自拍亚洲精品老妇| 国产精品久久久久久av不卡| 久久亚洲国产成人精品v| 免费av毛片视频| 99热全是精品| 精品熟女少妇av免费看| av视频在线观看入口| 夜夜爽夜夜爽视频| 国产精品av视频在线免费观看| 成人性生交大片免费视频hd| 国产高清不卡午夜福利| 日本免费a在线| 成人美女网站在线观看视频| 久久久久精品久久久久真实原创| 秋霞在线观看毛片| 国产成人精品久久久久久| 成人漫画全彩无遮挡| 亚洲国产精品国产精品| 国产真实伦视频高清在线观看| 毛片女人毛片| 一个人看视频在线观看www免费| 亚洲av免费在线观看| 啦啦啦啦在线视频资源| 国产午夜福利久久久久久| 嫩草影院精品99| 九九爱精品视频在线观看| 人妻夜夜爽99麻豆av| 亚洲最大成人av| 高清在线视频一区二区三区 | 久久久国产成人免费| 国产伦一二天堂av在线观看| 精品国产三级普通话版| 久久草成人影院| 国产免费福利视频在线观看| 免费观看人在逋| 午夜久久久久精精品| 国产av在哪里看| 欧美xxxx黑人xx丫x性爽| 午夜免费男女啪啪视频观看| 国产精品,欧美在线| 久久久久性生活片| 一级黄色大片毛片| 久久精品久久精品一区二区三区| 又爽又黄a免费视频| 又粗又爽又猛毛片免费看| 国产精品久久久久久av不卡| 欧美激情国产日韩精品一区| 欧美一区二区国产精品久久精品| 99热这里只有是精品50| 亚洲av男天堂| 欧美极品一区二区三区四区| 2022亚洲国产成人精品| 国产在视频线精品| av天堂中文字幕网| 岛国在线免费视频观看| 欧美日韩国产亚洲二区| 亚洲精品国产成人久久av| 精品人妻一区二区三区麻豆| 国产一区二区在线观看日韩| 狂野欧美白嫩少妇大欣赏| 美女脱内裤让男人舔精品视频| 日本免费a在线| 国产精品野战在线观看| 国产成年人精品一区二区| 成人亚洲精品av一区二区| 日韩av不卡免费在线播放| 少妇猛男粗大的猛烈进出视频 | 日韩视频在线欧美| 国产成人免费观看mmmm| 欧美3d第一页| 亚洲精品日韩在线中文字幕| 久久久久网色| 亚洲中文字幕日韩| 国产精品久久视频播放| 女人久久www免费人成看片 | 身体一侧抽搐| 大香蕉97超碰在线| 美女脱内裤让男人舔精品视频| 久久99热这里只频精品6学生 | 国产精品一区二区三区四区久久| 亚洲真实伦在线观看| 大香蕉97超碰在线| 夜夜看夜夜爽夜夜摸| 久久6这里有精品| 国产又色又爽无遮挡免| 高清av免费在线| 国产在视频线在精品| 久久国产乱子免费精品| 久久久国产成人精品二区| 久久久亚洲精品成人影院| 插阴视频在线观看视频| 欧美日韩综合久久久久久| 我要搜黄色片| 尾随美女入室| 欧美性猛交黑人性爽| 天堂网av新在线| 亚洲美女视频黄频| 中文字幕制服av| 亚洲精品aⅴ在线观看| 免费观看人在逋| 禁无遮挡网站| av在线天堂中文字幕| 国产精品久久久久久精品电影小说 | 男女国产视频网站| 少妇人妻精品综合一区二区| 日日摸夜夜添夜夜添av毛片| 日韩中字成人| 美女内射精品一级片tv| av在线老鸭窝| 91久久精品电影网| 国产黄色小视频在线观看| 一级黄片播放器| 欧美+日韩+精品| 黄片wwwwww| 一级毛片我不卡| 日本一本二区三区精品| 中文天堂在线官网| a级毛片免费高清观看在线播放| 国产精品一二三区在线看| 欧美人与善性xxx| 久久国产乱子免费精品| 卡戴珊不雅视频在线播放| av线在线观看网站| 久久久久久九九精品二区国产| 一级二级三级毛片免费看| 成人鲁丝片一二三区免费| 天天一区二区日本电影三级| 日本色播在线视频| 在线观看美女被高潮喷水网站| 亚洲久久久久久中文字幕| 天天一区二区日本电影三级| 热99在线观看视频| 精品一区二区三区视频在线| 国产伦在线观看视频一区| 最新中文字幕久久久久| 国产亚洲精品久久久com| 午夜视频国产福利| 国产精品福利在线免费观看| 69av精品久久久久久| 又黄又爽又刺激的免费视频.| 国产精品久久久久久久电影| АⅤ资源中文在线天堂| 亚洲精品乱码久久久v下载方式| 亚洲第一区二区三区不卡| 久久久久久久午夜电影| 免费看av在线观看网站| 少妇的逼水好多| 国产极品天堂在线| 搡老妇女老女人老熟妇| 中文字幕制服av| 国产成人精品久久久久久| 国产一区二区亚洲精品在线观看| 久99久视频精品免费| 国产免费男女视频| 精品一区二区免费观看| 日韩 亚洲 欧美在线| 国产毛片a区久久久久| 天堂中文最新版在线下载 | 日本免费a在线| 91精品国产九色| 全区人妻精品视频| 一夜夜www| 男女那种视频在线观看| 亚洲av中文av极速乱| 成人漫画全彩无遮挡| 毛片女人毛片| 九九爱精品视频在线观看| 神马国产精品三级电影在线观看| av卡一久久| 亚洲人成网站在线观看播放| 精品一区二区免费观看| 中文字幕精品亚洲无线码一区| 夫妻性生交免费视频一级片| 丰满乱子伦码专区| av线在线观看网站| 国产真实伦视频高清在线观看| 人妻少妇偷人精品九色| www.av在线官网国产| 精品久久久久久久久av| 色尼玛亚洲综合影院| 国产免费一级a男人的天堂| 亚洲最大成人av| 欧美极品一区二区三区四区| 夜夜看夜夜爽夜夜摸| 成人漫画全彩无遮挡| 国产精品麻豆人妻色哟哟久久 | 人人妻人人澡人人爽人人夜夜 | 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 男的添女的下面高潮视频| 久久国产乱子免费精品| 欧美xxxx黑人xx丫x性爽| 99久久精品一区二区三区| 九九久久精品国产亚洲av麻豆| 国产片特级美女逼逼视频| 精品久久久噜噜| 欧美日韩在线观看h| 男插女下体视频免费在线播放| 欧美日本视频| 日韩欧美三级三区| 亚洲av免费高清在线观看| 欧美最新免费一区二区三区| 一夜夜www| 亚洲国产日韩欧美精品在线观看| 国产一区二区在线av高清观看| 亚洲一区高清亚洲精品| 熟女电影av网| 少妇猛男粗大的猛烈进出视频 | 一边摸一边抽搐一进一小说| 内射极品少妇av片p| 久久99热这里只有精品18| 如何舔出高潮| 欧美一区二区精品小视频在线| 丝袜美腿在线中文| 国产成人freesex在线| 亚洲最大成人av| 国产精品电影一区二区三区| 亚洲美女视频黄频| 国产精品一二三区在线看| 亚洲精品aⅴ在线观看| 最后的刺客免费高清国语| 中文字幕av在线有码专区| 国产精品,欧美在线| 国产视频内射| 亚洲中文字幕日韩| 91在线精品国自产拍蜜月| 日韩在线高清观看一区二区三区| 久久99蜜桃精品久久| 搡女人真爽免费视频火全软件| 国产黄色视频一区二区在线观看 | 美女国产视频在线观看| 亚洲一区高清亚洲精品| 亚洲精品aⅴ在线观看| 卡戴珊不雅视频在线播放| 可以在线观看毛片的网站| 婷婷六月久久综合丁香| 又爽又黄无遮挡网站| 一区二区三区高清视频在线| 在现免费观看毛片| 欧美三级亚洲精品| 亚洲精品影视一区二区三区av| 久久久久免费精品人妻一区二区| 淫秽高清视频在线观看| 国产精品电影一区二区三区| 国产高清视频在线观看网站| 级片在线观看| 婷婷六月久久综合丁香| 一卡2卡三卡四卡精品乱码亚洲| 亚洲精品,欧美精品| 美女内射精品一级片tv| 国产伦在线观看视频一区| or卡值多少钱| a级一级毛片免费在线观看| 国产亚洲91精品色在线| 国产色婷婷99| 国产精品国产三级国产av玫瑰| av黄色大香蕉| 久久久成人免费电影| 啦啦啦啦在线视频资源| 久久人妻av系列| 欧美bdsm另类| 51国产日韩欧美| 美女脱内裤让男人舔精品视频| 色播亚洲综合网| 亚洲在久久综合| 青春草亚洲视频在线观看| 老女人水多毛片| 精品久久久久久电影网 | 白带黄色成豆腐渣| 一夜夜www| 特级一级黄色大片| 久久这里只有精品中国| 成人亚洲精品av一区二区| 国产一区二区三区av在线| 免费在线观看成人毛片| 亚洲中文字幕日韩| 99国产精品一区二区蜜桃av| 日韩欧美在线乱码| 久久久精品94久久精品| 国产一级毛片在线| 亚洲成人精品中文字幕电影| 青春草国产在线视频| 中文欧美无线码| 超碰97精品在线观看| 国产av在哪里看| 国产精品国产三级国产专区5o | 最近的中文字幕免费完整| 国产av码专区亚洲av| 男女视频在线观看网站免费| 麻豆一二三区av精品| 国产成人精品婷婷| 亚洲三级黄色毛片| 日本猛色少妇xxxxx猛交久久| 天堂网av新在线| 国产成人免费观看mmmm| 国产久久久一区二区三区| 免费看a级黄色片| 国产视频内射| 97人妻精品一区二区三区麻豆| 国产乱来视频区| 91午夜精品亚洲一区二区三区| 日韩 亚洲 欧美在线| 身体一侧抽搐| 草草在线视频免费看| 美女高潮的动态| 午夜福利在线在线| 国产美女午夜福利| 日本免费一区二区三区高清不卡| 1024手机看黄色片| 国产免费福利视频在线观看| 插阴视频在线观看视频| 亚洲色图av天堂| 久久韩国三级中文字幕| 日本熟妇午夜| 91久久精品国产一区二区成人| 亚洲国产精品合色在线| 嫩草影院新地址| 国产探花极品一区二区| 久久6这里有精品| 久久久久国产网址| 天美传媒精品一区二区| 两性午夜刺激爽爽歪歪视频在线观看| 精品免费久久久久久久清纯| 九草在线视频观看| 久久6这里有精品| 不卡视频在线观看欧美| or卡值多少钱| 中文资源天堂在线| 免费观看人在逋| 亚洲欧美日韩东京热| 日日干狠狠操夜夜爽| 国产精品一区二区三区四区免费观看| 日日干狠狠操夜夜爽| 国内精品一区二区在线观看| 69av精品久久久久久| 在线a可以看的网站| 久久久久九九精品影院| 亚洲av日韩在线播放| 国产真实乱freesex| 国产精品综合久久久久久久免费| 午夜激情欧美在线| 亚洲真实伦在线观看| 国产精品伦人一区二区| 26uuu在线亚洲综合色| 内地一区二区视频在线| 搡女人真爽免费视频火全软件| 在线播放国产精品三级| 日韩大片免费观看网站 | 成人毛片60女人毛片免费| 亚洲欧美日韩高清专用| 欧美一区二区精品小视频在线| 少妇人妻精品综合一区二区| 免费看a级黄色片| 欧美成人精品欧美一级黄| 热99在线观看视频| 中文在线观看免费www的网站| 国产一区亚洲一区在线观看| 大又大粗又爽又黄少妇毛片口| 九九热线精品视视频播放| 18+在线观看网站| 综合色av麻豆| 日韩成人伦理影院| 超碰av人人做人人爽久久| 两性午夜刺激爽爽歪歪视频在线观看| 女人十人毛片免费观看3o分钟| 色网站视频免费| 久久久久网色| 婷婷色麻豆天堂久久 | 色噜噜av男人的天堂激情| 别揉我奶头 嗯啊视频| 亚州av有码| 麻豆一二三区av精品| 国产精品女同一区二区软件| 国内精品一区二区在线观看| 成人无遮挡网站| 国产高潮美女av| 99热这里只有是精品在线观看| 婷婷色综合大香蕉| 亚洲图色成人| 国产一区亚洲一区在线观看| 麻豆一二三区av精品| 日本免费a在线| 黄色日韩在线| 99久久精品国产国产毛片| 亚洲av中文av极速乱| 女的被弄到高潮叫床怎么办| 婷婷六月久久综合丁香| 国产一区有黄有色的免费视频 | 久久久久精品久久久久真实原创| 丝袜喷水一区| 综合色av麻豆| 亚洲在久久综合| 午夜老司机福利剧场| 亚洲av日韩在线播放| 三级国产精品片| 久久这里只有精品中国| 亚洲在线观看片| 亚洲最大成人av| 欧美不卡视频在线免费观看| 亚洲综合精品二区| 天天躁夜夜躁狠狠久久av| 国产视频内射| 久久这里只有精品中国| 成人综合一区亚洲| 国产精品女同一区二区软件| 国产精品一区二区在线观看99 | av视频在线观看入口| 夜夜看夜夜爽夜夜摸| 非洲黑人性xxxx精品又粗又长| 少妇被粗大猛烈的视频| 国产日韩欧美在线精品| 床上黄色一级片| 久久久久久久久中文| 国产黄片视频在线免费观看| 免费无遮挡裸体视频| 在现免费观看毛片| 亚洲成人精品中文字幕电影| 内地一区二区视频在线| 高清av免费在线| 日本熟妇午夜| 禁无遮挡网站| 只有这里有精品99| 亚洲美女搞黄在线观看|