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

    Reconstruction and interpretation of photon Doppler velocimetry spectrum for ejecta particles from shock-loaded sample in vacuum?

    2021-06-26 03:04:12XiaoFengShi石曉峰DongJunMa馬東軍SonglinDang黨松琳ZongQiangMa馬宗強(qiáng)HaiQuanSun孫海權(quán)AnMinHe何安民andPeiWang王裴
    Chinese Physics B 2021年6期
    關(guān)鍵詞:安民海權(quán)馬東

    Xiao-Feng Shi(石曉峰) Dong-Jun Ma(馬東軍) Song-lin Dang(黨松琳) Zong-Qiang Ma(馬宗強(qiáng))Hai-Quan Sun(孫海權(quán)) An-Min He(何安民) and Pei Wang(王裴)

    1Institute of Applied Physical and Computational Mathematics,Beijing 100094,China

    2Jiangxi University of Applied Science,Nanchang 330103,China

    3Center for Applied Physics and Technology,Peking University,Beijing 100871,China

    Keywords: ejecta,photon Doppler velocimetry,Monte–Carlo algorithm,light scattering

    1. Introduction

    The strong shock wave released from the metal–vacuum/gas interface may eject a great number of metal particles.[1–5]Most of these particles are of micrometer-scale in size. This phenomenon of ejecta,or microjetting,was first observed by Kormeret al.in a plane impact experiment in the 1950s.[5]Later, the researchers in Los Alamos National Laboratory studied the production of particle ejection using the PHERMEX radiographic facility.[6]In recent decades, extensive investigations on particle ejection have been performed because of its important role in many scientific and engineering fields,including explosion damage,[7]pyrotechnics,[8]and inertial confinement fusion.[9,10]Many experimental approaches have attempted to measure the ejection production,such as the Asay foil,[3,11]foam recovery,[12]piezoelectric probes,[4,13]Fraunhofer holography,[14,15]x-ray/proton radiography,[2,16]Mie scattering,[16,17]and photon Doppler velocimetry (PDV).[18–24]In these approaches, Asay foil and piezoelectric probe are mainly used to measure the particles’momentum distribution, x-ray and proton radiography describe the space density of ejecta, foam recovery collects the total ejected mass, and Fraunhofer holography and Mie scattering measure the particles’diameter. These approaches can only measure some of these ejecta parameters. To reveal the full particle field of ejecta,multiple measurement approaches must be equipped. However, in real-world conditions, these approaches are hard to apply simultaneously because of limitations on the experimental space or configuration. Recently,PDV has attracted considerable attention[20,21,24]owing to its ability to recover the total area mass and the distributions of particle velocity and diameter at the same time. In addition,the light path of PDV is rather concise and its application is convenient. In some complex experimental configurations,PDV may be the only approach that can measure the ejecta particles.

    A standard PDV setup is shown in Fig. 1. The photodetector records a mixture of reference and backscattering light. The reference wave is in the carrier frequency, and the backscattering wave from ejecta particles has a shifted frequency due to the Doppler effect. The interference of the two light waves in the photodetector leads to temporal beats of light intensity. The beat signal consists of a large number of harmonics with different amplitudes and phases. The heterodyne signal may change according to variations in the particles’position and velocity. A discrete Fourier transform is applied to sweep the beats over time, giving a two-dimensional spectrogram on the “frequency/velocity–time” plane. In the spectrogram, the brightness of each point represents the corresponding spectral amplitude. The spectrogram is composed of the integral of all particles’scattering effects. Hence,interpreting the spectrogram in detail remains a challenging task.

    Fig. 1. Standard PDV setup: 1 laser; 2 reference light; 3 incident light; 4 ejecta;5 metal plate;6 shock;7 detonation;8 backscattering light;9 optical circulator;10 photodetector;11 photoelectric signal;12 PDV spectrogram;13 instantaneous spectrum.

    There have been several studies on the interpretation of the PDV spectrum. Buttler[25]used the spectrogram boundaries to determine the velocities of the spike and bubble of Richtmyer–Meshkov instability in loaded metal surface. The evolution of the PDV spectrogram in a gas environment was discussed by Sunet al.,[22]and the upper boundary of the spectrogram was used to obtain the particle size by considering aerodynamic deceleration effects. Fedorovet al.[23]discussed the influence of different particle sizes on the spectrogram boundary in further detail. Recently, Franzkowiaket al.[21]and Andriyashet al.[20,24]reconstructed the light field of ejecta and obtained the simulated PDV spectrum using single- and multiple-scattering theory,respectively. They varied the particles’parameters and fitted the simulated PDV spectrum to the experimental data. In this way, the particle velocity profile,diameter,and total area mass were recovered. Andriyashet al.considered the aerodynamic deceleration effects in a gas environment,whereas Franzkowiaket al. only discussed the case of a vacuum.

    Franzkowaiket al.[21]and Andriyashet al.[20,24]proposed similar approaches for recovering the ejecta parameters from the PDV spectrum through reconstruction and then fitting. However, some assumptions were introduced in the reconstruction of the light field.Franzkowaiket al.assumed that only backscattering light was present, while Andriyashet al.set the light scattering direction to be uniform and random in space. These assumptions affect the accuracy of the spectrum reconstruction, and thus influence the recovery of the ejecta parameters. The fitting model is another factor that affects the interpretation of the PDV spectrum. Different convergence criteria may produce different results. The quantitative relationship between the ejecta parameters and the PDV spectrum remains unclear. Hence,it is difficult to obtain definite ejecta parameters from the PDV spectrum. These issues provide the motivation for the present work.

    In this study, we improve the reconstruction method of the ejecta light field, and propose a novel model for extracting the ejecta parameters directly from the PDV spectrum.Mie theory, which gives a rigorous mathematical solution to Maxwell’s equations,is applied to calculate the light scattering effects,and a Monte–Carlo(MC)algorithm is used to describe the light transport process realistically. This reconstruction method provides a high-fidelity simulation for the PDV spectrum.The procedure is discussed in detail in Section 2.The influence of the ejecta parameters on the PDV parameters is then explored through MC simulations in Subsection 3.1. In Subsection 3.2, we propose an optical model that reveals the relationships between the PDV spectrum characteristics and the ejecta parameters. With this model,the ejecta parameters can be extracted directly from the PDV spectrum,instead of fitting to experimental data. In Subsection 3.3, the estimated ejecta parameters from an experimental PDV spectrum are verified against those measured by a piezoelectric probe. Finally, the conclusions to this study are presented in Section 4.

    2. Reconstruction of PDV spectrum

    2.1. Theoretical background

    The photodetector records reference and backscattering light waves. The scattering process of incident light is illustrated in Fig.2.The scattering light from the ejecta is governed by the superposition of waves propagating in the ejection particles along different light pathsi:

    whereEbsandEiare the electric vectors of total and partial scattering waves,respectively.

    The light intensity measured by the detector can be represented as

    whereE2randE2idenote the intensity of the reference and scattering light signals,respectively. The third term represents the heterodyne beats between the reference and scattering light,and the last term represents the heterodyne beats between the different scattering light paths. The Fourier transform ofI(t)is determined by the relation

    where|E|is the amplitude of the light wave field. In Eq.(3),only the third term appearing in Eq. (2) remains. This is because the frequencies of the first two terms are too high to be measured by the detector and the value of the last term is much smaller than that of the third term.ωris the carrier frequency,andωiis the Doppler-shifted frequency,which corresponds to a sequence of scattering events along pathi:

    wherecis the speed of light,nk,iandnsk,iare the directions of wave propagation before and after scattering by particlek,andvk,iis the velocity of this particle.

    Fig.2. Multiple scattering of light waves in ejection particles.

    To reconstruct the PDV spectrogram[Eq.(3)],the key is to obtain|Ei|andωi,i.e.,the detailed scattering process in the particles. For multiple-particle systems, the scattering field can be described by the transport equation[26,27]

    whereIis the light intensity in the field, which depends on both the detection positionrand the directionn.σsandκare the coefficients of scattering and absorption, respectively.p(n,n')is the scattering phase function. The right-hand side of this equation represents the contributions of scattering light from other positions.

    At the boundaries of the ejection, the light intensity has the form

    wheren0is the direction of incident light, which is usually perpendicular to the free surface.I0andIbsare the intensities of incident light and backscattering light, respectively, which correspond to the input and output of the transport equation.

    2.2. Monte–Carlo algorithm

    Andriyashet al.[20,24]used the discrete ordinate method to solve the transport equation[Eq.(5)]. In this paper,a more convenient and accurate method of the MC algorithm is applied to calculate the scattering effects. Compared with the discrete ordinate method, MC algorithm transform the complex light transport problem into multiple photons scattering problem. The detail scattering process can be considered in MC algorithm.

    In the MC algorithm,the incident light is assumed to be a great number of photons.When passing through random granular media,only part of the photons can penetrate.The proportion of permeable photons is approximated by Beer–Lambert’s law[17]

    whereLis the thickness of the medium andτdenotes the inverse extinction length,given by

    whereNis the number of particles per unit volume,Kextis the light extinction coefficient(determined by the particle diameter,light wavelength,and metal relative refraction index),is the mean cross-section area of the particles,diis the diameter of particlei,andAsis the light exposure area.

    For particles in the light exposure area,the total mass has the form

    wherem0is the total area mass of ejection andρ0is the density of the particle material.

    Combining Eqs.(8)and(9),we can rewrite Eq.(7)as

    The photons staying in the medium are scattered or absorbed by particles. The probabilities of scattering and absorption are calculated by the formula

    where the scattering coefficientKscaand the absorption coefficientKabsare calculated by Mie theory[28]as

    whereαis a dimensionless particle diameter parameter,α=πd/λ,λis the light wavelength, andan,bnare Mie coefficients. It is clear thatps+pa=1.

    If the photon is absorbed by particles, it will completely disappear and be converted into the particle’s internal energy.If the photon is scattered, its propagation direction and frequency will change,as shown in Fig.3. The phase function of the scattering angleθis calculated by the formula

    where Pnand P(1)nare Legendre and first-order associative Legendre functions,respectively.

    Fig.3. Light scattering on a particle.

    The incident light is assumed to be non-polarized,so the azimuth angle?after scattering obeys the uniform random distribution

    After scattering, the scattering angle and azimuth angle are added to the original light direction. The new direction cosine=[x,y,z]has the form

    whereu=[ux,uy,uz] is the original direction cosine. When|uz|≈1,the direction cosineis calculated by the formula

    After scattering, the frequency of the photon will have changed. The new frequency of scattering light has the form

    wherevis the velocity of the last particle that the photon left,is the current particle velocity,nis the photon direction,andωis the photon frequency before scattering.

    After passing through the entire ejection layer, few photons reach the free surface. An ideal diffuse reflection is assumed for these photons. After reflection,the space anglesθand?are uniformly random in [π/2,π] and [0,2π], respectively.

    Because the ejection is blocked by the free surface,eventually all the photons are either absorbed by particles or backscattered out from the top of the ejection. The frequency shifts of these“out”photons are summarized as the theoretical PDV spectrogram:

    wheren0is the number of initial photons andnoutis the distribution of out photons in terms of their frequency.

    The detailed steps of the calculation procedure are as follows:

    (i) First,the initial conditions of the photons and particles are set,such as the number and frequency of photons,and the diameter, velocity, and position of the particles. The photons start at the top of the ejection and then move towards the free surface.

    (ii) The step sizes of all photons are set to the same and equal to one thousandth of the height of the ejecta.

    (iii) In one iteration, all photons take one step in the direction of their propagation. Some photons may penetrate the current ejection layer, and the proportionpris determined by Eq. (10). For each photon, a random number is generated in(0,1). If the random number is less thenpr,the corresponding photon travels over the ejection layer boundary successfully.Otherwise,the corresponding photon is absorbed or scattered by particles in the ejection layer.

    (iv) For the photons that remain in the ejection layer,we use Eq. (11) to determine whether they are scattered or absorbed. If the photon is scattered, the direction change is calculated by Eqs. (13)–(17). Because the phase function of the scattering angle is very complex,an acceptance–rejection method is applied.The new frequency of the scattered photons is determined by Eq.(18).

    (v) Overall,if the photon travels across the ejection layer boundary,its position is updated;if the photon is scattered,its direction,frequency,and position are updated;if the photon is absorbed, it is labeled as such and removed from subsequent calculations.

    (vi) After updating the state of the photons, we check which of them have reached the free surface or left through the top of the ejection. For all photons that have reached the free surface,the ideal diffuse reflection is applied. If any photons have left the ejection, they are labeled accordingly and removed from subsequent calculations.

    (vii) Steps (iii)–(vi) are repeated until all photons have been absorbed or have left the ejection. The frequency shifts of outgoing photons are summarized as the spectrogram.

    2.3. Particle models

    The MC algorithm indicates that the PDV spectrum is related to the particle sized, velocityv, positionz, and numberN(i.e., total area massm0). This algorithm can be applied in cases where these parameters are completely random. In real situations, however, the particles of the ejecta usually satisfy certain distributions in terms of velocity and diameter.[2,14,15,17,29–32]For the sake of discussion, these assumptions are applied in this paper. Previous studies[2,29]indicate that the initial velocities of particles in the ejecta can be approximated by an exponential law

    wherevfsis the velocity of the free surface andβis the velocity distribution coefficient. Under this exponential law, most of the particles are located in the low-velocity region,which is near the free surface.βdetermines the non-uniformity of this distribution.

    In this paper,we only consider the ejecta in a vacuum environment. After being ejected, the particles retain an almost constant velocity and the ejecta expands in a self-similar manner over time. The particle positionzis only related to its initial velocityvand ejection timete,z=vte. The corresponding PDV spectrum exhibits slight changes over time.[21,33]

    The particle size distribution is assumed to obey a lognormal law[15,31]

    whereσis the width of the distribution anddmis the median diameter. These parameters depend on the roughness of the metal surface, shock-induced breakout pressure, and surrounding gas properties. Whenσ=0,the function becomes a Dirac equation and all of the particles have the same diameter.Obviously,this is the ideal situation. The particle distribution can also be described by a power law,[14,30,32]but this description may be invalid in the range of small particle sizes (less than 10 μm).[15]In this paper, the particle size is assumed to be independent of its velocity.

    With these assumptions, the determining factors of the PDV spectrum change to the velocity profile coefficientβ,total area massm0, median diameterdm, and size distribution widthσ. The aim of this paper is to discuss the influence of these parameters on the PDV spectrum, and to explore how they can be extracted from the PDV spectrum most accurately.

    2.4. Convergence and comparison

    The accuracy of the MC algorithm mainly depends on the initial number of photons. Theoretical PDV spectra with 104,105, 106, and 107initial photons are shown in Fig. 4. The calculation assumes a vacuum environment and the particles distribution assumptions are applied. In this case, the PDV spectra have a single peak. As the initial number of photons increases,the spectrum curves tend to be smooth. The difference between the spectra with 106and 107initial photons is very slight. Thus, 107initial photons are applied in the following calculations.

    Fig. 4. Theoretical PDV spectra with different initial numbers of photons.The calculation assumes the ejection of Sn particles in a vacuum environment. The particle velocities obey an exponential distribution(β =10)and the particle sizes follow a log-normal distribution(dm=1.5μm,σ =0.5).The total area mass is 20 mg/cm2. ωfs is the Doppler frequency shift corresponding to the free surface velocity, ωfs =2ωr·vfs/c. The probing wavelength is λ =1550 nm.

    The high initial number of photons leads to considerable computational cost. For the case of 107photons,a single-core CPU requires approximately 3 h to determine the spectrum.GPUs can be applied to accelerate the calculation. Although the frequency of GPU processors is much lower than that of CPUs,GPUs contain hundreds or thousands of stream processors that can work simultaneously. The acceleration ratio of a GPU compared to a CPU is shown in Fig.5. The Intel Xeon W-2102 CPU (frequency 2.9 GHz) and two GPUs (Quadro P600 and Nvidia GTX960) are applied. As the initial number of photons increases,the acceleration ratio of the GPUs is enhanced. For the case of 107photons, the acceleration ratio reaches a factor of 8 for the Quadro P600 and a factor of 20 for the Nvidia GTX960. Because there are many judgment events in the procedure,and the GPUs have few logical units,it is difficult to improve the acceleration ratio with these GPUs.Thus,the Nvidia GTX960, which requires approximately 500 s to compute each case,is used in the following calculations.

    Fig.5. Acceleration ratio of GPU calculation compared with CPU for different initial numbers of photons. The CPU is an Intel Xeon W-2102 and its basic frequency is 2.9 GHz. The Quadro P600 GPU has 384 stream processors; the clock speed of each processor is about 1.3 GHz. The Nvidia GTX960 GPU has 1024 stream processors;the clock speed of each processor is about 1.1 GHz.

    The PDV spectra simulated by the present procedure are compared with those reported by Andriyashet al. and Franzkowiaket al.in Fig.6. We use the equivalent ejecta area mass and particle size instead of the transport optical thickness used by Andriyashet al.With the uniform scattering assumption,our simulation(case 2)is almost the same as that of Andriyashet al.(case 3),which validates the adequacy of the present numerical method. However, when the Mie scattering theory is applied,there is a remarkable difference between the present procedure(case 1)and the results of Andriyashet al. (case 3) and Franzkowiaket al. (case 4). The difference with Andriyashet al.is mainly in the low-velocity part. This is because the change in the scattering phase function has a great influence on the multiple scattering, which is the main form of scattering in the low-velocity dense part. The difference with Franzkowiaket al. is in the location of the spectrum peak. Franzkowiaket al. applied the single-scattering theory and assumed that all of the light scattered backward.This implies that the optical thickness is overestimated,and so little light would reach the deep region of the ejecta. Thus,the spectrum moves toward high velocities. These differences in spectra indicate that the scattering assumption may introduce some reconstruction inaccuracy that cannot be neglected.

    Fig. 6. Comparison of PDV spectra calculated by different reconstruction methods. case 1: MC + Mie scattering theory (proposed procedure); case2: MC + uniform scattering assumption; case 3: Discrete coordinates +uniform scattering assumption(Andriyash et al.); case 4: Single scattering theory(Franzkowiak et al.).The data for case 3 were extracted directly from the paper of Andriyash et al.The calculations were carried out for a transport scattering thickness of τtr =10, which corresponds to m0 =10.7 mg/cm2,dm =1.5 μm, and σ =0.5. The material is Sn and the ejection velocity profile has β =8.

    3. Interpretation of PDV spectrum

    3.1. Influence of ejecta parameters

    The results of numerical calculations that demonstrate the sensitivity of the PDV spectrum to changes in the ejecta parameters (β,m0,dm, andσ) are presented in Figs. 7–10.The PDV spectra were simulated using the MC algorithm described in the previous section for Sn particles in a vacuum environment. The initial ejecta parameters were set toβ=10,m0=20 mg/cm2,dm=1.5μm,andσ=0.5. In each figure,one of the parameters changes and the others remain constant.

    Fig. 7. Simulated PDV spectra with different velocity coefficients. The calculation was carried out for Sn particles in a vacuum environment. The total area mass was 20 mg/cm2. The log-normal distribution(dm=1.5μm,σ =0.5)was applied to the particle sizes.

    The simulated PDV spectra with different values of the velocity coefficientβare shown in Fig. 7. With an increase inβ,the spectrum peak moves towards the low velocities and its magnitude increases. Furthermore,the spectrum shape becomes sharper and the high-velocity part of the spectrum becomes invisible. The coefficientβdetermines the distribution of particles in the ejecta. With largerβ, fewer particles are located at the top of the ejecta and the incident light can penetrate deeper. This results in the movement of the spectrum peak and a decrease in the observability of high-velocity particles.

    The simulated PDV spectra with different values of the total area massm0are shown in Fig. 8. The changes in the spectra can be divided into two sections. Whenm0≥

    10 mg/cm2, the spectrum displays a single peak. With a decrease inm0, this peak moves to the left, and its magnitude and slope exhibit slight changes. Whenm0≤5 mg/cm2, a new peak appears around the free surface, and the spectrum exhibits a double-peak shape. A smaller area mass produces a more remarkable new peak. Form0=2 mg/cm2, the original peak disappears and the spectrum again exhibits a single peak. The double-peak spectrum has been observed in previous experiments[18,34]and simulations,[24,35]and is the result of the direct exposure of incident light at the free surface. The detail explanation is shown in Subsection 3.2.

    Fig.8. Simulated PDV spectra with different total area mass. The area mass unit is mg/cm2. The calculation was carried out for Sn particles in a vacuum environment. The velocity coefficient β =10 and the size coefficients dm=1.5μm,σ =0.5.

    Fig. 9. Simulated PDV spectra with different particle median diameters.The diameter unit isμm. The calculation was carried out for Sn particles in a vacuum environment. The total area mass was 20 mg/cm2. The velocity coefficient β =10 mg/cm2 and the size coefficient σ =0.5.

    There are two parameters that determine the particle size distribution — the median diameterdmand the distribution widthσ. Their influence on the PDV spectrum is illustrated in Figs.9 and 10,respectively.dmandσexhibit similar effects:asdmorσincreases,the original peak of the spectrum moves towards the low velocities and the peak value decreases. A new peak then appears in the position of the free surface and the original peak gradually attenuates.This change in the form of the spectrum peak is similar to that for the area mass.

    Fig. 10. Simulated PDV spectra with different particle size coefficients σ.The calculation was carried out for Sn particles in a vacuum environment.The total area mass was 20 mg/cm2. The velocity coefficient β =10 and the median diameter dm=1.5μm.

    3.2. Theoretical optical model

    The simulations described above using the MC algorithm provide a qualitative understanding of the influence of the ejecta parameters on the PDV spectrum. However, how to solve the reverse problem,i.e., extracting the ejecta parameters from the PDV spectrum, remains unclear. To obtain the quantitative relationships between the ejecta parameters and the characteristics of the PDV spectrum,we introduce the single-scattering theory. In this theory,the light is assumed to be scattered only once, and the scattering direction is always backward. With this assumption,the light direction is always parallel to the motion of the particles. The frequency shift of the light is proportional to the particle velocity,ω=2ω0·v/c.The PDV spectrum can be expressed in terms of velocity,I().

    The above calculations have shown that the single scattering leads to an overestimation of the optical thickness. Here,we assume that the extinction process of particles only includes the backscattering and absorption effects, and the forward scattering is ignored. With this assumption, the PDV spectrum is calculated by the formula

    The exponent in Eq. (22) denotes the extinction effects,

    where the coefficient 2 signifies the back and forth of light in the ejection process. The integral represents the contributions from the extinction of particles above the layer of velocity:

    where the coefficientgbackis 0.5–0.7 for particle diameters of 1 μm–10 μm. Whengback=1, the present model reduces to that of Franzkowiaket al.

    This equation provides the theoretical form for determining the PDV spectrum from the ejecta parameters in a vacuum environment. In this formula, the PDV spectrum is proportional to the product of the velocity profile and the extinction term. These two parts are illustrated in Fig.11. As the velocity decreases, the corresponding ejecta position moves closer to the free surface,and the incident light becomes weaker because of particle extinction. However, the particles become dense deeper within the ejecta, and this enlarges the crosssection area of scattering. With the contribution of these two parts,the PDV spectrum exhibits a single peak shape.

    Fig.11. Extinction ratio and particle distribution with respect to velocity.

    Fig. 12. Simulated PDV spectrum by MC algorithm and single-scattering(SS)model.The area mass unit is mg/cm2.The particle settings are β =10,dm=1.5μm,σ =0.5. The backscattering coefficient is gback=0.67.

    The simulated PDV spectra of the present model are compared with those of the MC algorithm in Fig. 12. With the correction of the backscattering coefficient, there is a good agreement between the results, both in the main peak position and the curve shape. This agreement is resulted by the slight multiple scattering effects in ejecta field. According to the Mie theory, most of scattering light is forward with the particle diameter 1 μm–10 μm. So multiple scattering rarely occurs,and the single scattering theory can describe the PDV spectrum very well. However, in the case of a small ejecta mass (m0=5 mg/cm2), the present model cannot simulate the peak around the free surface. From the Eq.(22),the spectrum value in a certain velocity is related to the backscattering section and the local light intensity. In the free surface, most photons are scattered backward and the backscattering section is much larger than that of the particles,which results the new spectrum peak near the free surface. In the present model,the backscattering of incident light on the free surface is ignored,so there is no spectrum peak near the free surface in the case of a small ejecta mass(m0=5 mg/cm2). Even so,we mainly consider the information supplied by the original peak of the PDV spectrum in this paper. Thus,this defect has only a very slight influence on the accuracy of the present model.

    We now analyze the spectrum function[Eq.(28)]. First,we take its derivative

    In summary,the relationships between the ejecta parameters and the characteristics of the PDV spectrum have the form

    where the independent variable is normalized by the velocity of the free surface,vfs.

    It is worth nothing that the relationship (Eq. (35)) has some limitation. When the optical thickness is sufficiently small and the original peak disappears,it is hard to obtain the correct spectrum peak value and curvature. In this case, the ejecta parameter extraction method will be invalid.

    Fig.13. Peak positions of PDV spectrum calculated by MC algorithm and theoretical formula. Two optical thickness are considered,where the square points denote τ0 =23.4 and the circular points denote τ0 =8.6. The corresponding ejecta parameters are m0=20 mg/cm2,dm=1.5μm,σ =0.5 and m0=10 mg/cm2,dm=2.0μm,σ =0.5,respectively.

    Fig.14. Peak values of PDV spectrum calculated by MC algorithm and theoretical formula. Two optical thicknesses are considered,where the square points denote τ0 =23.4 and the circular points denote τ0 =8.6. The corresponding ejecta parameters are m0=20 mg/cm2,dm=1.5μm,σ =0.5 and m0=10 mg/cm2,dm=2.0μm,σ =0.5,respectively.

    Fig. 15. Relative curvature of PDV spectrum at the peak calculated by MC algorithm and theoretical formula. Two optical thickness are considered, where the square points denote τ0 =23.4 and the circular points denote τ0 =8.6. The corresponding ejecta parameters are m0 =20 mg/cm2,dm =1.5 μm, σ =0.5 and m0 =10 mg/cm2, dm =2.0 μm, σ =0.5, respectively.

    Figures 13–15 compare these theoretical relationships with the MC simulation results. Two optical thickness and a large range of velocity profiles are considered. In the cases ofτ0=23.4, the PDV spectra have single peak, and in the case ofτ0=8.6,the PDV spectra have double-peak. They are two typical types of PDV spectra. It can be observed that the theoretical relationships largely conform to the MC simulations in these cases,which verifies the present model to some extent.

    3.3. Experimental verification

    In the above relationships,the peak value of the spectrum is difficult to use in the analysis of PDV experiments. In the experiments,the PDV spectrum is scaled by the reference light intensity,probe reception,photoelectric conversion efficiency,and circuit amplification factor, among other factors. Additional PDV experiments are required to calibrate this scaled factor.For a single PDV vacuum experiment,only the velocity profileβand optical thicknessτ0of the ejecta can be extracted.If there is an additional particle granularity measurement,the ejecta area massm0can also be determined.

    Fig. 16. (a) PDV spectrogram extracted from ejecta experiment of Franzkowiak et al.[21] and(b)second derivative(deriv.) of the smooth data.

    A set of ejecta PDV experiments performed by Franzkowiaket al.[21]was used to verify the present theoretical model. The experiment was carried out in a vacuum environment using Sn material with the surface machined into 60 μm×8 μm grooves. The shock-induced breakout pressure wasPSB=28 GPa. The velocity of the free surface was found to be approximately 2013 m/s. We extracted the PDV spectrum from the experimental spectrogram over the period 0.2μs–0.8μs,as shown in Fig.16(a).The PDV data were then averaged and smoothed using the low-pass filtering of the fast Fourier transform. We converted the spectrum units[dBm]to volts and then took the second derivative to give the smoothed PDV spectrum shown in Fig. 16(b). The peak of this spectrum is located at/vfs=1.32 and the corresponding relative curvature is approximately?110. Combined with Eq. (35),this suggests a velocity profile coefficient ofβ=10.5 and an optical thickness ofτ0=28.79. Depending on the frequency of the filter, the errors of the extracted parameters are within 10%.

    Schaueret al.[31]conducted a Mie-scattering experiment with similar conditions, where the surface roughness was 50 μm×8 μm and the breakout pressure was about 30 GPa.The particle size distribution was measured to bedm=0.6μm,σ=0.5. Using this data, the total area mass was determined to bem0=7.5 mg/cm2.

    In their PDV experiment, Franzkowiaket al. simultaneously measured the area mass with respect to velocity using a piezoelectric probe. The PDV spectrum and area mass given by our estimations and their experiments are compared in Figs.17 and 18,respectively. These two results are in good agreement,which verifies the present theoretical model.

    Fig.17. Comparison of the experimental and simulated PDV spectra.

    Fig.18.Comparison of the area mass and velocity profile between the piezoelectric probe measurement[21] and our estimation.

    4. Conclusion

    This paper has discussed the PDV spectrum of ejecta particles from shock-loaded samples in a vacuum. A GPUaccelerated MC algorithm that rebuilds the PDV spectrum for the ejecta particles has been proposed,and Mie theory was applied to describe the scattering process. Compared with the reconstruction methods of Andriyashet al. and Franzkowiaket al., a reasonable scattering model is the key to simulating the PDV spectrum accurately. The influences of particle velocity profile,particle size and ejecta area mass on PDV spectrum are discussed,and two typical types of the PDV spectrum are found: single-peak and double-peak. For a quantitative analysis, a corrected single-scattering model is proposed for deriving the relationships between the ejecta parameters and the characteristics of the PDV spectrum. Based on the relationships,the ejecta parametersβ,d,andm0can be extracted from the PDV spectrum directly. The theoretical relationships and estimated parameters are found to be in good agreement with the MC simulations and PDV experiments of ejecta in a vacuum environment.

    The present theoretical model still has some limitation.The particles considered in this paper are in the diameter of 1 μm–10 μm, and the multiple scattering effects are very slight. If the particle diameter is smaller than this range, the multiple scattering effects must be considered. The spectrum reconstruction algorithm is still available, but the simplified theoretical model may be valid. In addition, the parameter extraction method is based on the original main peak. If the optical thickness is sufficiently small and the original peak disappears,the relationships between ejecta parameters and spectrum information will be hard to be utilized.How to determine the particle size in the PDV experiment is another unsolved issue. In future work,more experimental results will be used to verify the present model, and the PDV spectrum in gas environment will be discussed in an attempt to recover more comprehensive quantities of the ejecta.

    猜你喜歡
    安民海權(quán)馬東
    捐 款
    躬耕(2023年1期)2023-03-07 01:03:25
    THE EXISTENCE AND NON-EXISTENCE OFSIGN-CHANGING SOLUTIONS TO BI-HARMONIC EQUATIONS WITH A p-LAPLACIAN*
    Multiple solutions and hysteresis in the flows driven by surface with antisymmetric velocity profile?
    易安民聲
    易安民聲
    馬東生 作品
    晚清政府的海權(quán)意識(shí)與海軍實(shí)踐
    甲午戰(zhàn)爭與中國海權(quán)
    三坊七巷之安民巷
    馬漢及其『海權(quán)論』
    軍事歷史(1993年6期)1993-08-16 02:18:44
    手机成人av网站| 久久精品91无色码中文字幕| 国产高清视频在线播放一区| 可以在线观看毛片的网站| 99热6这里只有精品| 亚洲aⅴ乱码一区二区在线播放| 国产精品香港三级国产av潘金莲| 婷婷亚洲欧美| 91在线观看av| 亚洲男人的天堂狠狠| av在线天堂中文字幕| 丁香欧美五月| 动漫黄色视频在线观看| 天堂影院成人在线观看| 日本撒尿小便嘘嘘汇集6| 久久九九热精品免费| 69av精品久久久久久| 国产亚洲欧美在线一区二区| 亚洲18禁久久av| 国产99白浆流出| 99热精品在线国产| 午夜精品在线福利| 国产精品爽爽va在线观看网站| 亚洲国产精品sss在线观看| 国产精品99久久久久久久久| 级片在线观看| 欧美日韩福利视频一区二区| 国产精品日韩av在线免费观看| 国产午夜精品久久久久久一区二区三区 | 亚洲av成人精品一区久久| 黄色片一级片一级黄色片| 69人妻影院| 老汉色av国产亚洲站长工具| 亚洲欧美日韩无卡精品| 老司机在亚洲福利影院| 精品国内亚洲2022精品成人| av黄色大香蕉| 最新美女视频免费是黄的| 麻豆成人午夜福利视频| 99国产综合亚洲精品| 国产欧美日韩一区二区三| 国内少妇人妻偷人精品xxx网站| 成人永久免费在线观看视频| 一卡2卡三卡四卡精品乱码亚洲| 嫁个100分男人电影在线观看| 噜噜噜噜噜久久久久久91| 男女下面进入的视频免费午夜| 国产一区二区三区视频了| 亚洲无线观看免费| 精品久久久久久久久久久久久| 真人一进一出gif抽搐免费| 亚洲av一区综合| 国产欧美日韩一区二区三| 亚洲 国产 在线| 一个人看的www免费观看视频| 美女被艹到高潮喷水动态| 亚洲av电影不卡..在线观看| 校园春色视频在线观看| 琪琪午夜伦伦电影理论片6080| 国产一区二区在线观看日韩 | 波多野结衣高清无吗| 757午夜福利合集在线观看| 欧美日本视频| 国产av麻豆久久久久久久| 亚洲在线观看片| 国产亚洲欧美98| 给我免费播放毛片高清在线观看| 日韩欧美免费精品| 日本精品一区二区三区蜜桃| 中文字幕人成人乱码亚洲影| 久久久久久九九精品二区国产| 黄色成人免费大全| 99久久九九国产精品国产免费| 久久久久久久久中文| 啦啦啦韩国在线观看视频| 成人欧美大片| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 国产一区在线观看成人免费| www.色视频.com| 一a级毛片在线观看| 夜夜看夜夜爽夜夜摸| 国产精品乱码一区二三区的特点| 人妻丰满熟妇av一区二区三区| 亚洲精品456在线播放app | 国产精品98久久久久久宅男小说| 免费观看精品视频网站| 亚洲av不卡在线观看| 久久草成人影院| 欧美黄色淫秽网站| 99国产精品一区二区蜜桃av| 国产一区二区亚洲精品在线观看| 国内揄拍国产精品人妻在线| 神马国产精品三级电影在线观看| av天堂中文字幕网| 亚洲一区二区三区不卡视频| 啦啦啦观看免费观看视频高清| 亚洲avbb在线观看| 日本一本二区三区精品| 淫妇啪啪啪对白视频| 日韩欧美国产在线观看| 午夜激情欧美在线| 日韩高清综合在线| 天美传媒精品一区二区| 国产av不卡久久| 亚洲激情在线av| 在线视频色国产色| 日韩欧美国产一区二区入口| 亚洲乱码一区二区免费版| av中文乱码字幕在线| 中文字幕人成人乱码亚洲影| 久久精品国产自在天天线| 女人高潮潮喷娇喘18禁视频| 我的老师免费观看完整版| 又紧又爽又黄一区二区| 熟妇人妻久久中文字幕3abv| 18+在线观看网站| 窝窝影院91人妻| 黑人欧美特级aaaaaa片| 内地一区二区视频在线| 亚洲国产中文字幕在线视频| 午夜久久久久精精品| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 色综合亚洲欧美另类图片| 亚洲精品日韩av片在线观看 | 亚洲精品粉嫩美女一区| 舔av片在线| 免费人成视频x8x8入口观看| 日韩亚洲欧美综合| 国产亚洲av嫩草精品影院| 我要搜黄色片| 中文字幕熟女人妻在线| 国产欧美日韩精品亚洲av| 国产欧美日韩一区二区三| 午夜a级毛片| 丁香欧美五月| 日本五十路高清| 国产极品精品免费视频能看的| 精品不卡国产一区二区三区| 黄色日韩在线| 久久久久免费精品人妻一区二区| 韩国av一区二区三区四区| 香蕉久久夜色| 波野结衣二区三区在线 | 一区二区三区高清视频在线| 在线观看免费视频日本深夜| av福利片在线观看| 久久国产精品影院| 国产精品98久久久久久宅男小说| 成年版毛片免费区| 又黄又爽又免费观看的视频| 白带黄色成豆腐渣| 成年版毛片免费区| 色哟哟哟哟哟哟| 又紧又爽又黄一区二区| 黄色日韩在线| 亚洲精品一卡2卡三卡4卡5卡| 激情在线观看视频在线高清| 国产伦一二天堂av在线观看| 一区二区三区国产精品乱码| 国产欧美日韩精品一区二区| 国产色婷婷99| 99热精品在线国产| 精品国内亚洲2022精品成人| 亚洲欧美日韩高清在线视频| 99精品在免费线老司机午夜| 99久久精品一区二区三区| 内地一区二区视频在线| 国产精品久久电影中文字幕| 欧美成人免费av一区二区三区| 久久精品国产99精品国产亚洲性色| 国产视频内射| 欧美黄色淫秽网站| 少妇的丰满在线观看| 国产成人欧美在线观看| 欧美丝袜亚洲另类 | 动漫黄色视频在线观看| 亚洲美女黄片视频| 亚洲成人中文字幕在线播放| 搡女人真爽免费视频火全软件 | 欧美区成人在线视频| a级一级毛片免费在线观看| 蜜桃亚洲精品一区二区三区| 热99re8久久精品国产| 成人亚洲精品av一区二区| 亚洲久久久久久中文字幕| 丝袜美腿在线中文| 一级作爱视频免费观看| 最近最新中文字幕大全电影3| 好看av亚洲va欧美ⅴa在| 白带黄色成豆腐渣| 波野结衣二区三区在线 | 国内精品美女久久久久久| 欧美一区二区精品小视频在线| 国产精品久久久久久亚洲av鲁大| 99久久精品热视频| 俄罗斯特黄特色一大片| 亚洲国产色片| 亚洲男人的天堂狠狠| 日本一本二区三区精品| 麻豆一二三区av精品| 亚洲国产欧美网| 手机成人av网站| 两个人的视频大全免费| av天堂在线播放| 国产私拍福利视频在线观看| 天天添夜夜摸| 亚洲aⅴ乱码一区二区在线播放| 一边摸一边抽搐一进一小说| 成人av在线播放网站| xxxwww97欧美| 真实男女啪啪啪动态图| 国产亚洲精品av在线| 成人特级黄色片久久久久久久| 国产精品久久久久久亚洲av鲁大| 久久久久亚洲av毛片大全| 99国产极品粉嫩在线观看| 精品国产三级普通话版| 丰满乱子伦码专区| 69av精品久久久久久| 久久久色成人| 99久国产av精品| 真人做人爱边吃奶动态| 香蕉久久夜色| h日本视频在线播放| 全区人妻精品视频| 中文字幕av成人在线电影| 村上凉子中文字幕在线| avwww免费| 天堂动漫精品| 乱人视频在线观看| 国产精品久久久久久人妻精品电影| 99riav亚洲国产免费| 成人性生交大片免费视频hd| 桃色一区二区三区在线观看| 欧美另类亚洲清纯唯美| 欧美大码av| 看片在线看免费视频| 日本熟妇午夜| 亚洲国产高清在线一区二区三| 日韩人妻高清精品专区| 男女午夜视频在线观看| 亚洲美女视频黄频| 免费高清视频大片| 我要搜黄色片| 一本久久中文字幕| 国产中年淑女户外野战色| 亚洲精品456在线播放app | 欧美又色又爽又黄视频| 日韩欧美精品v在线| 丰满的人妻完整版| 久久久国产精品麻豆| 最近最新免费中文字幕在线| 熟女人妻精品中文字幕| 亚洲精品456在线播放app | 久久久久久久久久黄片| 日日干狠狠操夜夜爽| 成人av在线播放网站| 中出人妻视频一区二区| 国产激情偷乱视频一区二区| 日本成人三级电影网站| 欧美日本亚洲视频在线播放| 国产av不卡久久| 欧美zozozo另类| 国产淫片久久久久久久久 | 亚洲熟妇中文字幕五十中出| 免费看美女性在线毛片视频| 99riav亚洲国产免费| 国内精品久久久久久久电影| 久久伊人香网站| 看片在线看免费视频| 国产乱人伦免费视频| 久久久久久九九精品二区国产| 一a级毛片在线观看| 亚洲熟妇熟女久久| 欧美成人性av电影在线观看| 国产欧美日韩一区二区精品| 精品日产1卡2卡| 99热精品在线国产| 又黄又爽又免费观看的视频| 亚洲成人久久性| 久久精品国产自在天天线| 99在线人妻在线中文字幕| 69av精品久久久久久| 99精品久久久久人妻精品| 欧美极品一区二区三区四区| 狂野欧美激情性xxxx| 亚洲第一欧美日韩一区二区三区| 国产亚洲欧美98| 国产v大片淫在线免费观看| www.熟女人妻精品国产| 美女高潮喷水抽搐中文字幕| 91在线精品国自产拍蜜月 | 久久国产精品影院| 午夜福利欧美成人| 久久久久久久精品吃奶| 真实男女啪啪啪动态图| 12—13女人毛片做爰片一| 精品免费久久久久久久清纯| 黄色丝袜av网址大全| 欧美区成人在线视频| 欧美bdsm另类| 色精品久久人妻99蜜桃| 国产精品 欧美亚洲| 精品不卡国产一区二区三区| 久久久久久久久大av| 女人十人毛片免费观看3o分钟| 在线看三级毛片| 一区二区三区国产精品乱码| av在线蜜桃| 国产精品一区二区三区四区久久| 欧美乱妇无乱码| 老司机福利观看| 国产精品亚洲美女久久久| 每晚都被弄得嗷嗷叫到高潮| 久久久久久久久大av| 国产亚洲精品久久久com| 亚洲不卡免费看| 国产成+人综合+亚洲专区| 最好的美女福利视频网| 亚洲狠狠婷婷综合久久图片| 1024手机看黄色片| svipshipincom国产片| 白带黄色成豆腐渣| 一进一出好大好爽视频| 欧美成人一区二区免费高清观看| 真实男女啪啪啪动态图| 白带黄色成豆腐渣| 国产成人系列免费观看| 国产一区二区激情短视频| 成人无遮挡网站| 成熟少妇高潮喷水视频| 看黄色毛片网站| 我的老师免费观看完整版| 欧美午夜高清在线| 69av精品久久久久久| 最新在线观看一区二区三区| 美女被艹到高潮喷水动态| 久久精品91蜜桃| 久久久成人免费电影| 一夜夜www| 国产视频内射| 一级毛片高清免费大全| 老汉色∧v一级毛片| 18禁国产床啪视频网站| 国产精品久久久久久人妻精品电影| 成人鲁丝片一二三区免费| 午夜免费男女啪啪视频观看 | 白带黄色成豆腐渣| 亚洲自拍偷在线| 成年女人永久免费观看视频| 在线观看免费午夜福利视频| av专区在线播放| av天堂在线播放| 少妇裸体淫交视频免费看高清| 欧美日韩福利视频一区二区| 国内精品久久久久精免费| 久久国产精品影院| 99久久精品热视频| 欧美成人a在线观看| 欧美中文综合在线视频| 亚洲av五月六月丁香网| 男人舔女人下体高潮全视频| 久久久久久久精品吃奶| 欧美中文综合在线视频| 又黄又粗又硬又大视频| 亚洲国产精品sss在线观看| 免费搜索国产男女视频| 国产久久久一区二区三区| 小说图片视频综合网站| 男人舔女人下体高潮全视频| 亚洲精品影视一区二区三区av| 丰满人妻熟妇乱又伦精品不卡| 午夜两性在线视频| 亚洲精品日韩av片在线观看 | 51国产日韩欧美| 香蕉久久夜色| 久久久久久久久大av| 国产欧美日韩一区二区精品| 制服人妻中文乱码| 国产视频一区二区在线看| 男女做爰动态图高潮gif福利片| 免费搜索国产男女视频| 精品久久久久久成人av| 亚洲国产精品合色在线| avwww免费| 久久精品91无色码中文字幕| xxxwww97欧美| 国产免费男女视频| 亚洲成人中文字幕在线播放| 狠狠狠狠99中文字幕| 亚洲av不卡在线观看| 白带黄色成豆腐渣| 黄色女人牲交| 嫁个100分男人电影在线观看| 最近最新中文字幕大全电影3| 精品国产亚洲在线| 久久精品影院6| bbb黄色大片| 欧美高清成人免费视频www| 亚洲精品一区av在线观看| 国内毛片毛片毛片毛片毛片| 国产v大片淫在线免费观看| 亚洲精品乱码久久久v下载方式 | 69av精品久久久久久| 欧美日韩亚洲国产一区二区在线观看| 国内精品一区二区在线观看| 男人舔女人下体高潮全视频| 久久天躁狠狠躁夜夜2o2o| 看片在线看免费视频| 日韩 欧美 亚洲 中文字幕| 欧美黑人巨大hd| 国产高清videossex| 国产免费一级a男人的天堂| 免费大片18禁| 精品人妻一区二区三区麻豆 | 亚洲欧美精品综合久久99| 国语自产精品视频在线第100页| 午夜免费观看网址| 欧美色视频一区免费| 国产亚洲精品综合一区在线观看| 91av网一区二区| 欧美成人a在线观看| 久久久精品大字幕| 成人国产一区最新在线观看| 我要搜黄色片| 色尼玛亚洲综合影院| 国产乱人视频| 亚洲黑人精品在线| av专区在线播放| 国产一区在线观看成人免费| 国产亚洲精品久久久com| 看片在线看免费视频| 美女免费视频网站| 日韩欧美精品免费久久 | 丰满人妻熟妇乱又伦精品不卡| 亚洲电影在线观看av| 国产一区二区在线av高清观看| e午夜精品久久久久久久| 国产精品国产高清国产av| 91麻豆精品激情在线观看国产| 日韩高清综合在线| 国产激情偷乱视频一区二区| 午夜影院日韩av| 色精品久久人妻99蜜桃| 一级毛片女人18水好多| 国产av在哪里看| 欧美成人a在线观看| 一级毛片女人18水好多| 真人一进一出gif抽搐免费| 国产欧美日韩精品一区二区| 国产午夜精品久久久久久一区二区三区 | 真人做人爱边吃奶动态| 亚洲自拍偷在线| www.www免费av| 看片在线看免费视频| 国产不卡一卡二| 可以在线观看的亚洲视频| 黄色成人免费大全| 久久久精品欧美日韩精品| 国产在线精品亚洲第一网站| 97人妻精品一区二区三区麻豆| 亚洲av免费高清在线观看| 久久久久国产精品人妻aⅴ院| 久久久久久久午夜电影| 国产成人欧美在线观看| 国产免费一级a男人的天堂| 色播亚洲综合网| 少妇的逼水好多| 精品国产三级普通话版| 99精品久久久久人妻精品| 国产精品国产高清国产av| 有码 亚洲区| 尤物成人国产欧美一区二区三区| 婷婷精品国产亚洲av在线| 18美女黄网站色大片免费观看| 欧美高清成人免费视频www| 亚洲激情在线av| 成年免费大片在线观看| 久久精品91蜜桃| 免费搜索国产男女视频| 成人精品一区二区免费| 国产蜜桃级精品一区二区三区| 老司机午夜福利在线观看视频| 精品国产亚洲在线| 精品免费久久久久久久清纯| 国产欧美日韩精品一区二区| 久久99热这里只有精品18| av片东京热男人的天堂| 精品一区二区三区人妻视频| 一夜夜www| 成人一区二区视频在线观看| 午夜a级毛片| 最新在线观看一区二区三区| 久久精品91蜜桃| 国产高清有码在线观看视频| 真人做人爱边吃奶动态| 欧美午夜高清在线| 婷婷丁香在线五月| 日韩精品青青久久久久久| 美女cb高潮喷水在线观看| 久久精品国产亚洲av涩爱 | 国产精品久久久久久久久免 | 窝窝影院91人妻| 欧美午夜高清在线| 中出人妻视频一区二区| 日韩欧美国产一区二区入口| 男人的好看免费观看在线视频| 又黄又粗又硬又大视频| 国产国拍精品亚洲av在线观看 | 国产亚洲精品av在线| 亚洲,欧美精品.| 狂野欧美激情性xxxx| 久久久久久人人人人人| 嫩草影院入口| 欧美高清成人免费视频www| 久久中文看片网| 久久精品国产99精品国产亚洲性色| 不卡一级毛片| 欧美黑人巨大hd| 国产亚洲精品一区二区www| 免费高清视频大片| 国产久久久一区二区三区| 女人高潮潮喷娇喘18禁视频| 一本一本综合久久| 国语自产精品视频在线第100页| 久久久精品大字幕| 亚洲熟妇中文字幕五十中出| 亚洲人成网站在线播| 午夜免费成人在线视频| 在线观看日韩欧美| 亚洲av成人不卡在线观看播放网| 欧美日韩福利视频一区二区| 在线观看av片永久免费下载| 成人性生交大片免费视频hd| 免费av不卡在线播放| 国内毛片毛片毛片毛片毛片| 啦啦啦韩国在线观看视频| 一进一出好大好爽视频| 亚洲第一欧美日韩一区二区三区| 亚洲 国产 在线| 国产精品嫩草影院av在线观看 | 久久国产乱子伦精品免费另类| 99在线视频只有这里精品首页| av福利片在线观看| 亚洲七黄色美女视频| 国产亚洲精品综合一区在线观看| 天堂动漫精品| 午夜久久久久精精品| 在线观看一区二区三区| 国产69精品久久久久777片| 日韩人妻高清精品专区| 亚洲av免费高清在线观看| 欧美最新免费一区二区三区 | 国产淫片久久久久久久久 | 高潮久久久久久久久久久不卡| 久久久成人免费电影| 在线看三级毛片| 国产高清视频在线播放一区| 久久性视频一级片| 国产老妇女一区| 天天躁日日操中文字幕| 亚洲av二区三区四区| 女人十人毛片免费观看3o分钟| 国产免费男女视频| 日韩欧美免费精品| 精品久久久久久久人妻蜜臀av| 变态另类丝袜制服| 亚洲av成人精品一区久久| 免费av毛片视频| 蜜桃久久精品国产亚洲av| 日韩 欧美 亚洲 中文字幕| 首页视频小说图片口味搜索| 神马国产精品三级电影在线观看| 麻豆成人午夜福利视频| 亚洲中文字幕一区二区三区有码在线看| 欧美成人一区二区免费高清观看| 真实男女啪啪啪动态图| 高清在线国产一区| 男女下面进入的视频免费午夜| 国产精品久久久久久久久免 | 久久久久精品国产欧美久久久| 一区福利在线观看| 国产又黄又爽又无遮挡在线| 五月玫瑰六月丁香| 成人欧美大片| 久久香蕉国产精品| 看黄色毛片网站| 久久久精品欧美日韩精品| 久久精品国产亚洲av涩爱 | 全区人妻精品视频| 午夜福利视频1000在线观看| 中文亚洲av片在线观看爽| 嫁个100分男人电影在线观看| 全区人妻精品视频| 国产成人影院久久av| 夜夜看夜夜爽夜夜摸| 国产亚洲av嫩草精品影院| 亚洲成人免费电影在线观看| 99riav亚洲国产免费| 精品国产亚洲在线| 欧美色视频一区免费| 久久久久久人人人人人| 国产麻豆成人av免费视频| 欧美又色又爽又黄视频| 免费在线观看亚洲国产| 亚洲美女黄片视频| 特大巨黑吊av在线直播| av欧美777| 中文字幕人妻熟人妻熟丝袜美 | 法律面前人人平等表现在哪些方面| 久久精品国产亚洲av香蕉五月| 欧美大码av| 亚洲国产欧美人成| 午夜免费男女啪啪视频观看 | 亚洲精品久久国产高清桃花|