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

    Axial acoustic radiation force on an elastic spherical shell near an impedance boundary for zero-order quasi-Bessel–Gauss beam?

    2021-05-06 08:55:00YuChenZang臧雨宸WeiJunLin林偉軍ChangSu蘇暢andPengFeiWu吳鵬飛
    Chinese Physics B 2021年4期
    關鍵詞:鵬飛

    Yu-Chen Zang(臧雨宸), Wei-Jun Lin(林偉軍),?, Chang Su(蘇暢), and Peng-Fei Wu(吳鵬飛)

    1Institute of Acoustics,Chinese Academy of Sciences,Beijing 100190,China

    2University of Chinese Academy of Sciences,Beijing 100049,China

    Keywords: acoustic radiation force,zero-order quasi-Bessel–Gauss beam,impedance boundary,elastic spherical shell

    1. Introduction

    The acoustic radiation force, arising from the momentum transfer between sound waves and scattering objects,has been a hot topic due to its potential applications in acoustic tweezers,[1–5]acoustic levitation,[6,7]and particle manipulation.[8–10]Since King[11]first set forth the concept of acoustic radiation force on particles and calculated the force experienced by a rigid sphere in travelling and standing wave field, there have been many relevant theoretical and experimental studies till now. A large number of these studies are concerned with the radiation force acting on spherical objects as well as the determination of the acoustic intensity through measuring the radiation force. Yosioka et al.[12]first developed a theory for the acoustic radiation force acting on a compressible sphere. Hasegawa et al.[13]successfully calculated the radiation force on an elastic sphere, showing that both the compressive and shear wave velocities play critical roles in determining the radiation force. Shortly after, they presented an analytical expression of acoustic radiation force on a sphere in a quasi-stationary field.[14]With the development of beamforming and ultrasonic transducer technology, Wu &Du[15]first considered the radiation force in a focused Gauss beam field, laying the foundation of developing single beam acoustic tweezers. A single beam acoustic tweezer for threedimensional(3D)trapping in a vortex beam was also presented by Baresch et al.,[16]which can selectively pick and control the position of a particle in three dimensions as well as achieve the controlled rotation of absorbing particles. Marston[17–20]and Mitri[21–28]systematically investigated the acoustic radiation force acting on a sphere exerted by both zero-order and high order Bessel beams.In recent years,considering the commonness of the spheroidal objects in actual practices such as the suspended droplets in air,some researchers have extended their scope of study to the spheroids in Bessel beams.[29–32]

    Based on all of these results, there is a growing interest among researchers to consider the acoustic radiation force on small particles near an impedance boundary instead of an unbounded field. Actually, the particle manipulation in practice are often performed in a bounded field especially in many biomedical applications. Miri et al.[33]added a vessel porous wall to the single spherical contrast agent shell and analyzed the effect of the impedance boundary on acoustic radiation force in detail. Wang et al.[34]theoretically and numerically studied the acoustic radiation force exerted by a standing wave on a rigid cylinder near a flat wall. Qiao et al.[35,36]developed a model to work out the acoustic radiation force acting on a fluid cylinder near an impedance boundary using the image theory. More recently,the same authors extended the research to the elastic cylinder case.[37]Analyses of the acoustic radiation force on a fluid sphere between two impedance boundaries were also performed by our group.[38]

    It is well known that micro-shells have been utilized by many researchers for drug and gene delivery systems.[39–41]Due to its non-immunogenity and noninvasiveness, the ultrasound-mediated microbubble system has drawn wide attention from researchers in both the diagnostic and therapeutic applications of ultrasound contrast agents. To safely and accurately guide the contrast agent micro-shells towards the desired targeted positions,it is important to make clear of the mechanism of acoustic radiation force acting on this kind of structure. Hasegawa et al.[42]solved the acoustic radiation force exerted on both cylindrical and spherical shells in plane progressive waves. Later, Mitri[43,44]extended this theory to the case of elastic and viscoelastic spherical shells in plane standing waves, in which the sound absorption effect inside the elastic material was systematically studied. Based on the results in earlier works,our group also presented a theoretical model to compute the radiation force on a multilayered elastic shell in standing Bessel beams.[45]

    As a typical type of non-diffracting beams,Bessel beams have both the property of narrow beam-width and large depth of field, thus they have many potential applications in ultrasonic imaging and other fields. Given an infinite aperture, a zero-order Bessel beam can travel to infinity without changing its profile and amplitude. However, the aperture sizes of the beams are always finite, therefore only within the depth of field is the fundamental Bessel beam viewed as non-diffracting. With the radial profile given by Bessel–Gauss function, Bessel–Gauss beams can be seen as the simplest case of the physically realizable finite aperture Bessel beams.In 2003, the axicon-based Bessel–Gauss resonator with concave output coupler was firstly developed in optics.[46]Later,a simple and compact source was created by researchers to directly produce a Bessel–Gauss beam.[47]Since then,the nondiffracting property of Bessel–Gauss beams has been used to generate low-order harmonics in gas.[48]As for its applications in acoustics, the properties of the fundamental and second harmonic components were studied by Ding et al.[49]Some approximate analytical descriptions were also investigated with a finite aperture.[50,51]In 2017,Chen et al.[52]first combined the study of acoustic radiation force with quasi-Bessel–Gauss beams, theoretically figuring out the radiation force on a sphere in a progressive and standing zero-order quasi-Bessel–Gauss beam. They found that the energy concentration of Gauss component can make up for the shortfall of the relatively small amplitude of acoustic radiation force in Bessel beams.However,to the best of our knowledge,no work has been done dealing with the acoustic radiation force near a boundary in this case.This study aims to extend Chen’s results to a bounded space,providing a closed-form expression of the acoustic radiation force on an elastic spherical shell near an impedance interface for zero-order quasi-Bessel–Gauss beam.This model can be taken as the simulation of a contrast agent besides the human blood vessel wall. Admittedly, a cluster of micro-shells are more likely to be used in drug-loaded ultrasound applications. However, the study on a single shell serves as a necessary interim step to deepen our understanding of the radiation force effect. Also, a proper estimation of the radiation force is of great significance to avoid the shell’s fragmentation and attachment to the blood vessel wall. Some numerical simulations are also carried out with particular emphasis on: (i) the effects of the impedance boundary; (ii) the effects of the shell’s thickness; (iii) the half-cone angle and beam waist of the incident beam;(iv)the fluid medium filling the hollow space; (v)the viscosity of the coated elastic layer.For purposes of this study,it is assumed that the amplitudes of both the sound pressure and particle velocity are so small that the nonlinear wave propagation effects can be neglected.

    2. Theory

    The geometry of this problem is shown in Fig.1.A spherical shell of outer radius a and inner radius b is positioned in a non-viscous fluid and at distance d from the impedance boundary. The center of the shell coincides with the origin of both the Cartesian coordinate system(x,y,z)and the spherical coordinate system (r, θ, φ). From the geometric relationship between two systems, we can conclude that x=rsinθ cosφ,y=rsinθ sinφ, z=rcosθ. Let a zero-order quasi-Bessel–Gauss beam be incident upon the shell along the+z direction,which is perpendicular to the impedance boundary. As a complex solution of the paraxial equation, the standard Bessel–Gauss beam has been studied by previous researchers in detail,which features the non-diffracting advantage of Bessel beams and the energy focusing property of Gauss beams.[53,54]

    2.1. Expansion of incident beam in terms of spherical functions

    Assume that the field distribution at the waist of a zeroorder quasi-Bessel–Gauss beam is given by

    where ω is the angular frequency of the incident beam,α and w0are two constants, r⊥=rsinθ is the radial distance from the axis of the beam,and J0denotes the zero-order cylindrical Bessel function. We can imagine it as a Gauss beam “modulated” in amplitude by the Bessel function, or a zero-order Bessel function weighted by the Gaussian aperture. The two parameters α and w0determine the effective width and aperture of a zero-order quasi-Bessel–Gauss beam, respectively.Obviously, the zero-order Bessel beam and the Gauss beam are special cases of the zero-order quasi-Bessel–Gauss beam.It will degenerate into a basic Gauss beam when α approaches zero and will degenerate into a basic zero-order Bessel beam when w0approaches infinity. Substituting Eq. (1) into the Fresnel propagation formula,we can obtain the mathematical expression of the beam field at distance z as follows:[47]

    where k=ω/c0is the wave number in the fluid medium, c0is the compressional wave speed in the surrounding fluid, zRis the Rayleigh range, Φ(z)=kz ?α2z/2k ?arctan(z/zR), is the axial phase of the beam. Moreover,the spot size w(z)and the wave front curvature R(z)are the same as those for a basic Gauss beam.

    Strictly speaking, a quasi-Bessel–Gauss beam does not possess the non-diffracting property. In other words,the radial distribution depends on the propagation distance, which can be easily seen from Eq.(2). However,it can be approximated to be non-diffracting in the region near the focal plane. Based on this premise,we express the local approximation for the velocity potential of a quasi-Bessel–Gauss beam near the focal region at z=0 as[52]

    where φ0is the amplitude of the velocity potential,kz=kcosβ and kr=ksinβ are the axial and radial wave numbers,β is the half-cone angle of the Bessel beam component and W0is the beam waist of the Gauss beam component.

    In order to enhance the intuitive understanding of the zero-order quasi-Bessel–Gauss beam,the radial distribution of velocity potential at z=0 with different half-cone angles and beam waists are shown in Fig.2 based on Eq.(3). As can be seen from the plots,increasing the half-cone angle β leads to the decrease of the main lobe width.Moreover,the beam waist W0mainly determines how fast the acoustic energy attenuates in the radial direction. In other words, a smaller beam waist corresponds to a quasi-Bessel–Gauss beam with a stronger focusing characteristic.

    According to the finite series expansion theory, the incident wave can be expanded in spherical coordinates with orthogonal properties and finite series properties as[47]

    where jnis the n-th-order spherical Bessel function of the first kind,and Pnis the n-th-order Legendre function. The expanding coefficients Cncan be calculated by using the special value method[55]involving expansion with spherical harmonic functions. Substituting n by 2l and 2l+1 for even and odd n,respectively,the coefficients Cncan be expressed as[52]

    2.2. Scattering of the incident beam by a spherical shell near the boundary

    The reflection wave of the incident beam from the impedance boundary can be expressed in the spherical coordinate system as

    where Rs=(Z ?ρ0c0)/(Z+ρ0c0)is the reflecting coefficient of the impedance boundary, in which Z is the surface normal acoustic impedance of the medium where the image spherical shell is immersed, ρ0is the mass density of the fluid medium. It is noted that Rs=1 corresponds to the absolutely rigid boundary condition,while Rs=0 corresponds to the unbounded condition.

    The scattering wave velocity potential of the real and image shell can be respectively expressed as

    According to the theory of vector analysis, the particle displacement in the coated elastic material can be expressed as the sum of the gradient of a scalar potential and the curl of a vector potential as

    where φ2and ψ satisfy the wave equations for the compressional and shear waves filling the elastic material,respectively as

    where ?2denotes the Laplacian operator, kc= ω/ccand ks=ω/csrepresent the compressional and shear wave numbers,respectively. Their wave velocities can be calculated by

    where λ and μ denote the Lame coefficients, ρ is the mass density of the coated elastic layer. Hence, φ2and ψ can be expressed as

    where yndenotes the n-th-order spherical Bessel function of the second kind,bn,cn,dn,enare the coefficients determined by certain boundary conditions.

    As for the sound wave in the ideal fluid filling the hollow space,its velocity potential can be written as

    where k3is the wave number in the inner fluid medium.

    Obviously, the boundary conditions in this problem are set at both r=a and r=b: (i) The acoustic pressure in the fluid must be equal to the normal component stress in the solid;(ii)The normal component of the displacement must be continuous in the interfaces;(iii)The tangential component of the stress must be equal to zero in the interfaces.

    It seems that we can directly work out the problem of acoustic scattering through the boundary conditions above.However,the velocity potentials are expressed in different coordinate systems,presenting a challenge to our computing process.Using the classical form of the translational addition theorem for bi-spherical coordinates,we can express the scattering velocity potential from the image shell in spherical coordinate system(r,θ,φ)as follows:

    where the detailed expression of the parameter Qmnis given in Appendix B,which can also be found in Ref.[56].

    Equation(15)can also be rewritten in another form as

    where

    Therefore,the normal velocity and the total sound pressure can be expressed as

    Incorporating Eqs.(13),(14),(17),and(18)into the boundary conditions mentioned above,we can obtain a linear system of equations. The detailed expressions of these equations are presented in Appendix C,thus only a simplified matrix form of them is given here as follows:

    2.3. Acoustic radiation force on a sphere near the boundary

    Considering the symmetrical property of this problem,only the axial component of the radiation force acting on the spherical shell needs to be computed. Projecting Eq. (20) in the direction of z axis,we can obtain

    Therefore,the axial radiation force can be expressed by

    where

    After some mathematical calculations,the radiation force can finally be written in a simplified form as

    where

    in which Re and Im represent the real and imaginary parts of a complex function,respectively.

    3. Numerical simulations and analyses

    In order to more clearly illustrate the general behavior of the acoustic radiation force on an elastic spherical shell near an impedance boundary, several numerical simulations are provided in this section. The physical constants of all the materials used in the following simulations are listed in Table 1.Hereafter, the surrounding ideal fluid is assumed to be water.In all of our calculations,the infinite series in Eq.(25)are truncated at|an/a0|~10?6,which can well balance the accuracy and calculation efficiency.

    Material Density Compressible Shear(kg/m3) velocity(m/s) velocity(m/s)Water 1000 1480 ...Air 1.23 344 ...Aluminum 2700 6420 3040 Polymethyl methacrylate 1190 2690 1340

    3.1. Effects of the impedance boundary on the radiation force

    The effect of the impedance boundary on the radiation force function curves for a single elastic spherical shell is discussed in this section. Figure 4 plots the radiation force functions on an elastic shell with b/a = 0.75 and d/a = 2 for a zero-order quasi-Bessel–Gauss beam in the case W0=3λ,β =π/3. The fluid medium within the hollow space is assumed to be water. In Figs. 4(a) and 4(b), the elastic shell is supposed to be made of polymethyl methacrylate(PMMA)and aluminum, respectively, which are two common materials in biomedical ultrasound. The reflecting coefficient of the impedance boundary is set as Rs=0.1,0.4,0.7,and 1.0. For a blood vessel,the reflecting coefficient is estimated to be in the range 0.05~0.15 according to the measurement for the acoustic impedance of blood vessel walls.[57]However, a variable reflecting coefficient ranging from 0 to 1 is considered in the following calculations to better visualize the effect of changing Rson the radiation force. Also,the case when Rs=0.1 is especially studied. The results for the radiation force function in free space are also provided in the same figure for ease of comparison. In our simulation,the increment step of ka is selected as 0.01 so as to avoid missing the very narrow resonant peaks. Inspection of the results reflects the huge influence of the interface added to the problem. Compared with the results in free space, the radiation force function curves exhibit an obvious oscillatory property. For the case of a PMMA shell,a series of sharp peaks appear in the curves, which may correspond to the resonant frequencies of the elastic vibration in the shell. At these peaks, nearly all the acoustic energy is stored within the shell without leaking into the fluid. However, for an aluminum shell, the high-amplitude peaks disappear. This is mainly due to the larger discontinuity of the acoustic impedance at the water-aluminum interface than the water-PMMA interface,thus less acoustic energy can be transmitted into the elastic material. As the reflecting coefficient Rsincreases,the amplitudes of the peaks will be significantly enhanced. However,the variation of Rsdoes not shift the positions of maxima and minima in the curves. In addition, it is interesting to observe that for both spherical shells, the radiation force can be negative at selected ka values, which is another distinctive property compared with the results in an unbounded fluid medium. This phenomenon means that the shell can be driven towards the sound source, resulting from the reflection wave by the impedance boundary. An explicit explanation of the physical mechanism for the negative radiation force is that the scattering in the forward hemisphere is relatively stronger than in the backward hemisphere, which has been proved in earlier studies on Bessel beams.[17,58]

    To better understand the features of the acoustic radiation force on an elastic spherical shell for zero-order Bessel–Gauss beams,the total scattering velocity potential in the far-field region can be expressed in the term of form function as

    For an absolutely rigid sphere, the form function can be expressed as:

    The theoretical discussion above is applied to the PMMA spherical shell in Fig.4(a). The form function moduli of the total, rigid background, and pure resonance back scattering are depicted in Figs. 5(a), 5(b), and 5(c), respectively. It is shown that the rigid background scattering plots oscillate around|f|=1 when ka>2 but do not have any sharp peaks.The oscillations in Fig.5(b) are caused by the interferences between the specular contribution and Franz wave. The difference between the total and rigid background scattering form functions is due to the resonance of different vibration modes,as is clearly shown in the pure resonance backscattering form function plots. It is easy to find that the resonance scattering peaks in Fig.5(c)coincide well with the radiation force function peaks in Fig.4(a), which proves that all of these sharp peaks result from the resonant vibration modes. These resonant vibrations are closely related to the acoustic property of the shell and the surrounding fluid medium(water).

    In real applications, the spherical shells are often at various distances from the reflecting interface. Therefore, it is necessary to consider the radiation force functions with different distance-to-radius ratios. With Rs=0.5 and b/a=0.75,the radiation force functions on a water-filled PMMA shell at different d/a values are depicted in Fig.6. Other parameters about the incident beam are set as β =π/3 and W0=3 λ. It is noticeable that the extrema in the figures,corresponding to the resonance vibration modes of the shell,can appear as not only the maxima but also the minima. This phenomenon can also be found in the simulated results for a spherical shell in free space in Ref.[44].As the spherical shell moves away from the boundary, which corresponds to the increase of d/a, the amplitudes of the resonant peaks will decrease obviously at high frequencies. These phenomena prove that the shell–boundary distance is an important factor influencing the radiation force peaks.

    To more accurately evaluate the effects of distance-toradius ratio on the radiation force,the Ypversus d/a curves for a water-filled PMMA shell at ka=0.42 in the case of different Rsare shown in Fig.7,which can be regarded as the incidence of a 10-MHz beam upon a micro-shell of outer radius 10μm.All of the other parameters are set exactly the same as in Fig.6.It is easy to find that the peaks and dips are periodically distributed along the d/a axis. As has been found in Ref. [35],when the boundary is absolutely rigid (Rs=1), the radiation force vanishes at pressure nodes(kd =π/2, 3π/2, 5π/2,...)or antinodes(kd=0,π,2π,...)and gets its maximum value at the intermediate locations(kd=π/4,3π/4,5π/4,...).Therefore, the period of the Ypcurves along the d/a axis can be calculated as π/ka=7.48 given that ka=0.42, which consists well with the simulated plots. With the growth of Rs,the amplitude of Ypincreases,consisting well with the conclusion from Fig.4.

    3.2. Effects of the shell’s thickness on the radiation force

    In biomedical ultrasound applications, the drug-loaded micro-shells are often developed with various thicknesses to meet different demands in practice. To evaluate the effect of changing the spherical shell’s relative thickness, the computed results for the radiation force functions on a waterfilled PMMA shell with different relative thicknesses are depicted in Fig.8. With Rs= 0.5, d/a = 2, β = π/3, and W0= 3λ, the two subgraphs (a) and (b) plot the results in the case 0.6 ≤b/a ≤0.9 with increment ?(b/a) = 0.1 and 0.92 ≤b/a ≤0.98 with increment ?(b/a) = 0.02, respectively. It is interesting to find in Fig.8(a) that when the coated elastic layer is relatively thick (b/a<0.9), highamplitude peaks emerge at resonant frequencies, which correspond to different ka values with the variation of b/a.However, the resonant peaks disappear in the simulated ka range when the value of b/a reaches 0.9. In fact, with the increase of b/a, the resonant vibration modes shift to high frequencies that are beyond our simulated range. Nevertheless,we are not going to study these peaks considering the severe attenuation of sound energy at high frequencies. The results for thinner shells in Fig.8(b)also exhibit no sharp resonant peaks,while a series of maxima and minima still exist.The increase of b/a leads to the decrease of Ypbut does not influence the positions of peaks and dips. It is also noteworthy that the attractive force can be obtained at appropriate ka and b/a values.

    3.3. Effects of the incident beam on the radiation force

    As has been mentioned above,a zero-order quasi-Bessel–Gauss beam is characterized by the half-cone angle β and the beam waist W0. Ypis simulated in the half-cone angle range of 0<β <π/2 with increment 0.01π, as well as in the ka range of 0

    It is also noteworthy that taking β =0 for the zero-order Bessel–Gauss beam yields to the results in the basic Gauss beam field. For the case of a PMMA shell, the value of Ypis likely to exceed 70 at selected ka values when β is relatively small. Admittedly, the high peaks in the simulated plots will definitely be damped once the viscosity is considered. However,a theoretical guide can still be provided for acoustic manipulation applications based on these peaks when a strong force is required to push or pull the target particles more easily. As β increases from 0 to π/2, the amplitude of Ypdecreases rapidly at most ka values. It can be expected that the force vanishes at β =π/2 based on the property of Legendre function that Pn(0)Pn+1(0)=0,which is perfectly confirmed by our simulated results.

    Next,we turn to evaluate how the variation of beam waist affects the radiation force. Unlike the half-cone angle β, the beam waist W0determines the Gauss component of a zeroorder quasi-Bessel–Gauss beam. Figure 10 shows the computed results for a PMMA shell and an aluminum shell in the range of 3λ

    3.4. Effects of the material within hollow space on the radiation force

    Gas-filled micro-shells are also very common in the applications of contrast agents, which play an important role in modern drug delivery and imaging systems. Assuming that the shell is filled with air instead of water, the simulated results for a PMMA shell and an aluminum shell with different reflecting coefficients are presented in Figs. 11(a) and 11(b),respectively. The parameters are set as b/a=0.75,W0=3λ,β = π/3, d = 2a, which are exactly the same as those in Fig.4 for ease of comparison. For the case of a PMMA shell, the number of sharp resonant peaks is reduced, but their amplitudes see a pronounced increase.This phenomenon reflects that the change of the interior core material can significantly affect the resonant frequencies. From the view point of manipulating target particles,an air-filled PMMA shell is a preferable choice compared with a water-filled PMMA shell.However, for an aluminum shell, there exists no obvious difference between the radiation force function curves for an airfilled shell and a water-filled shell. This is mainly because that only little energy can enter the interior core due to the serious mismatching of acoustic impedance between the adjacent layers. Therefore,the variation of the material filling the hollow space does not obviously affect the property of radiation force.

    3.5. Effects of the viscosity on the radiation force

    4. Conclusions

    Based on the image theory and the finite series method,a zero-order quasi-Bessel–Gauss beam has been expanded in terms of harmonic spherical functions. A closed-form expression of the acoustic radiation force on an elastic spherical shell near an impedance boundary for zero-order quasi-Bessel–Gauss beam is derived. On the basis of theoretical results, several corresponding simulated results are also provided. The oscillatory property and the negative radiation force at selected ka values are two main influences of the impedance boundary added to this problem. Some sharp resonant peaks appear in the results for a PMMA shell but not an aluminum shell. With the increase of Rs,high-amplitude radiation force can be obtained,while the positions of maxima and minima remain unchanged. When the shell moves along the beam axis,the radiation force function will vanish at pressure nodes and antinodes and increase to maximum at the intermediate locations. For a PMMA shell, a relatively thick coated layer is conducive to pushing or pulling particles through the strong radiation force. As the two parameters characterizing a zero-order quasi-Bessel–Gauss beam,the half-cone angle and the beam waist can both significantly influence the radiation force. Resulting from the energy focusing effect of the Gauss component, the radiation force function exhibits an obvious increase at small beam waists in the range ka>4. The viscosity of the elastic coated layer leads to the attenuation of sound energy and then the reduction of radiation force peaks at high ka values. In summary,the present work systematically investigated the acoustic radiation force on an elastic spherical shell near the boundary,which can be extended to the drug delivery and imaging in biomedical ultrasound.

    Appendix B:Expression ofQmn

    The paprameter Qmnis given by

    Appendix C:Expressions of Eq.(19)

    The detailed expressions of Eq.(19)are given as follows:

    猜你喜歡
    鵬飛
    樊應舉
    書香兩岸(2020年3期)2020-06-29 12:33:45
    漫畫
    Quality Control for Traditional Medicines - Chinese Crude Drugs
    為了避嫌
    雜文月刊(2019年18期)2019-12-04 08:30:40
    懲“前”毖“后”
    21世紀(2019年10期)2019-11-02 03:17:02
    執(zhí)“迷”不悟
    21世紀(2019年10期)2019-11-02 03:17:02
    舉賢
    21世紀(2019年9期)2019-10-12 06:33:44
    漫畫
    粗看“段”,細看“端”
    漫畫
    狠狠精品人妻久久久久久综合| 黑丝袜美女国产一区| 国产乱人偷精品视频| 精品少妇一区二区三区视频日本电影 | 搡老乐熟女国产| 亚洲伊人久久精品综合| 精品久久久久久电影网| 欧美 亚洲 国产 日韩一| 男女国产视频网站| 日韩中文字幕视频在线看片| 日韩精品有码人妻一区| 中文字幕人妻熟女乱码| 久久99蜜桃精品久久| 亚洲精品久久久久久婷婷小说| 亚洲精品自拍成人| 成人亚洲欧美一区二区av| 久久鲁丝午夜福利片| 巨乳人妻的诱惑在线观看| 成年av动漫网址| 免费黄色在线免费观看| 婷婷色综合www| 欧美亚洲 丝袜 人妻 在线| 激情视频va一区二区三区| 永久网站在线| 啦啦啦啦在线视频资源| 精品一区二区免费观看| 80岁老熟妇乱子伦牲交| 免费观看无遮挡的男女| 一本色道久久久久久精品综合| 91aial.com中文字幕在线观看| 制服丝袜香蕉在线| 岛国毛片在线播放| 久久婷婷青草| 丝袜喷水一区| www.自偷自拍.com| 黄色 视频免费看| 久久久久久免费高清国产稀缺| 亚洲精品成人av观看孕妇| 高清欧美精品videossex| 精品少妇久久久久久888优播| 精品久久蜜臀av无| 如日韩欧美国产精品一区二区三区| 国产成人精品久久二区二区91 | 亚洲国产精品成人久久小说| 亚洲 欧美一区二区三区| 成人影院久久| 国产激情久久老熟女| 亚洲av男天堂| 26uuu在线亚洲综合色| 久久免费观看电影| 大话2 男鬼变身卡| 美女脱内裤让男人舔精品视频| 欧美人与善性xxx| 青春草视频在线免费观看| 黄色 视频免费看| 免费看av在线观看网站| 美女国产高潮福利片在线看| 午夜日韩欧美国产| 国产精品一区二区在线不卡| 性少妇av在线| 国产成人精品无人区| 美女高潮到喷水免费观看| 免费看不卡的av| 天天影视国产精品| 九九爱精品视频在线观看| 黄色 视频免费看| 国产无遮挡羞羞视频在线观看| 午夜激情av网站| 一级毛片我不卡| 久久ye,这里只有精品| 一边亲一边摸免费视频| 一级毛片 在线播放| 日本vs欧美在线观看视频| 男人操女人黄网站| 日产精品乱码卡一卡2卡三| 高清黄色对白视频在线免费看| av.在线天堂| 日韩熟女老妇一区二区性免费视频| 久久精品aⅴ一区二区三区四区 | 女人久久www免费人成看片| 国产一区二区激情短视频 | 成年女人在线观看亚洲视频| 国产亚洲一区二区精品| 国产亚洲一区二区精品| 热re99久久精品国产66热6| 久久久精品区二区三区| 97精品久久久久久久久久精品| 男人爽女人下面视频在线观看| 99久国产av精品国产电影| 亚洲四区av| 一本大道久久a久久精品| 一区二区三区乱码不卡18| 国产成人av激情在线播放| 一级毛片电影观看| www.熟女人妻精品国产| 99热网站在线观看| 一本—道久久a久久精品蜜桃钙片| 中国三级夫妇交换| 女人精品久久久久毛片| 精品久久蜜臀av无| 又粗又硬又长又爽又黄的视频| 日韩三级伦理在线观看| 自拍欧美九色日韩亚洲蝌蚪91| 免费久久久久久久精品成人欧美视频| 久久久久精品久久久久真实原创| 国产精品无大码| 丝袜人妻中文字幕| 免费人妻精品一区二区三区视频| 香蕉国产在线看| 街头女战士在线观看网站| 男女下面插进去视频免费观看| 天堂中文最新版在线下载| 亚洲欧美精品自产自拍| 纵有疾风起免费观看全集完整版| 欧美日韩一区二区视频在线观看视频在线| 久久精品国产亚洲av高清一级| 男女啪啪激烈高潮av片| 伊人亚洲综合成人网| 制服诱惑二区| 一本—道久久a久久精品蜜桃钙片| 久久久久国产一级毛片高清牌| 日日撸夜夜添| 亚洲精品国产色婷婷电影| 精品少妇久久久久久888优播| 天堂8中文在线网| 亚洲国产最新在线播放| 2021少妇久久久久久久久久久| 99久久精品国产国产毛片| 国产 精品1| 老熟女久久久| 九九爱精品视频在线观看| 九九爱精品视频在线观看| 极品少妇高潮喷水抽搐| av在线app专区| 女人高潮潮喷娇喘18禁视频| 久久久a久久爽久久v久久| 国产野战对白在线观看| 亚洲少妇的诱惑av| 巨乳人妻的诱惑在线观看| 久久精品熟女亚洲av麻豆精品| av在线老鸭窝| 伦理电影大哥的女人| 看免费成人av毛片| 90打野战视频偷拍视频| 永久免费av网站大全| 国产极品天堂在线| 美女主播在线视频| 最新中文字幕久久久久| 菩萨蛮人人尽说江南好唐韦庄| 成人国产av品久久久| 日韩免费高清中文字幕av| 9191精品国产免费久久| 亚洲色图综合在线观看| av卡一久久| 美女国产高潮福利片在线看| 免费在线观看视频国产中文字幕亚洲 | 免费女性裸体啪啪无遮挡网站| 国产成人精品久久二区二区91 | videos熟女内射| av在线老鸭窝| 热99久久久久精品小说推荐| 老汉色av国产亚洲站长工具| 国产精品一区二区在线观看99| 熟女电影av网| 亚洲一区二区三区欧美精品| 美女大奶头黄色视频| 国产黄色免费在线视频| 欧美人与性动交α欧美精品济南到 | 最黄视频免费看| av网站在线播放免费| 午夜激情久久久久久久| 黑人猛操日本美女一级片| 国产一区二区 视频在线| 久久国产亚洲av麻豆专区| 日本欧美国产在线视频| 人妻系列 视频| 国产欧美日韩一区二区三区在线| 日韩制服丝袜自拍偷拍| 国产免费又黄又爽又色| 热99久久久久精品小说推荐| 欧美+日韩+精品| 久久精品人人爽人人爽视色| 母亲3免费完整高清在线观看 | 啦啦啦中文免费视频观看日本| 激情视频va一区二区三区| 免费黄网站久久成人精品| a级片在线免费高清观看视频| 亚洲少妇的诱惑av| 热99久久久久精品小说推荐| 99国产综合亚洲精品| 欧美日韩一级在线毛片| 黑人欧美特级aaaaaa片| 99热全是精品| 十八禁高潮呻吟视频| 久久av网站| 国产高清不卡午夜福利| av网站免费在线观看视频| av免费观看日本| 国产精品人妻久久久影院| 在线观看免费高清a一片| 亚洲欧美日韩另类电影网站| 丝袜喷水一区| 成人亚洲欧美一区二区av| 91国产中文字幕| 亚洲美女视频黄频| 中文精品一卡2卡3卡4更新| 久久青草综合色| 亚洲欧美精品综合一区二区三区 | 欧美成人午夜免费资源| 国产又爽黄色视频| 人妻 亚洲 视频| av有码第一页| av网站免费在线观看视频| 国产 一区精品| 国产成人91sexporn| 欧美少妇被猛烈插入视频| 久久国产亚洲av麻豆专区| 久久久久网色| 如何舔出高潮| 欧美激情 高清一区二区三区| 国产熟女欧美一区二区| 国产精品偷伦视频观看了| 女人精品久久久久毛片| 欧美日韩视频高清一区二区三区二| 好男人视频免费观看在线| 成年人午夜在线观看视频| 免费黄网站久久成人精品| 国产精品香港三级国产av潘金莲 | 国产熟女午夜一区二区三区| 又大又黄又爽视频免费| 免费av中文字幕在线| www日本在线高清视频| 美女中出高潮动态图| 亚洲欧美精品自产自拍| 色哟哟·www| 最近的中文字幕免费完整| 欧美亚洲日本最大视频资源| 国产一区有黄有色的免费视频| 欧美+日韩+精品| 午夜福利在线观看免费完整高清在| 日本-黄色视频高清免费观看| videosex国产| 丝袜美腿诱惑在线| 香蕉国产在线看| 寂寞人妻少妇视频99o| 赤兔流量卡办理| 一级爰片在线观看| 亚洲精品美女久久久久99蜜臀 | 国产av国产精品国产| 国产 精品1| 伊人久久国产一区二区| a级片在线免费高清观看视频| 狠狠婷婷综合久久久久久88av| 这个男人来自地球电影免费观看 | 99热国产这里只有精品6| √禁漫天堂资源中文www| 亚洲国产欧美网| 蜜桃国产av成人99| 欧美变态另类bdsm刘玥| 99国产综合亚洲精品| 男人操女人黄网站| 亚洲内射少妇av| 久久精品国产综合久久久| 亚洲五月色婷婷综合| 国产日韩欧美亚洲二区| 欧美人与善性xxx| 狠狠婷婷综合久久久久久88av| 如日韩欧美国产精品一区二区三区| 日韩大片免费观看网站| 人人妻人人添人人爽欧美一区卜| 热99国产精品久久久久久7| 国产爽快片一区二区三区| 精品酒店卫生间| 日韩精品有码人妻一区| 青青草视频在线视频观看| 三级国产精品片| 1024香蕉在线观看| 久久久久久久久久久久大奶| 免费不卡的大黄色大毛片视频在线观看| 国产一区二区 视频在线| 久久久久网色| 在线观看人妻少妇| 精品国产一区二区三区久久久樱花| 午夜av观看不卡| 久久午夜综合久久蜜桃| 欧美人与性动交α欧美软件| 久久精品国产a三级三级三级| 99国产精品免费福利视频| 国产精品 国内视频| 国产精品人妻久久久影院| 国产成人免费无遮挡视频| 亚洲欧美日韩另类电影网站| 久久久国产欧美日韩av| 捣出白浆h1v1| 亚洲国产精品一区二区三区在线| 国产精品成人在线| 免费av中文字幕在线| 国产片特级美女逼逼视频| 高清在线视频一区二区三区| 久久亚洲国产成人精品v| 中文字幕人妻丝袜制服| 日韩制服骚丝袜av| 免费观看性生交大片5| 一区二区三区激情视频| 人妻一区二区av| 纵有疾风起免费观看全集完整版| 日本vs欧美在线观看视频| 王馨瑶露胸无遮挡在线观看| 伊人久久大香线蕉亚洲五| 深夜精品福利| 日韩人妻精品一区2区三区| 精品午夜福利在线看| 亚洲av免费高清在线观看| 欧美97在线视频| a级毛片在线看网站| 男女边吃奶边做爰视频| 午夜91福利影院| 久久久久国产精品人妻一区二区| 少妇精品久久久久久久| 成人亚洲精品一区在线观看| 成人黄色视频免费在线看| 久久人人爽av亚洲精品天堂| 成人二区视频| kizo精华| 国产熟女午夜一区二区三区| 啦啦啦视频在线资源免费观看| 飞空精品影院首页| 久久狼人影院| 亚洲欧美一区二区三区黑人 | 亚洲欧美精品综合一区二区三区 | 美女脱内裤让男人舔精品视频| 欧美日韩一级在线毛片| 亚洲精品久久午夜乱码| 欧美日韩亚洲高清精品| 美女高潮到喷水免费观看| 国产精品成人在线| 国产男女内射视频| 春色校园在线视频观看| 亚洲av在线观看美女高潮| 久久久久国产网址| 男女下面插进去视频免费观看| 搡老乐熟女国产| 亚洲精品国产av蜜桃| 亚洲久久久国产精品| 国产黄频视频在线观看| 久久热在线av| 亚洲精品在线美女| 人人妻人人爽人人添夜夜欢视频| 男女边摸边吃奶| www.熟女人妻精品国产| 在线免费观看不下载黄p国产| 亚洲视频免费观看视频| 婷婷色麻豆天堂久久| 男女国产视频网站| 麻豆精品久久久久久蜜桃| 久久亚洲国产成人精品v| 国产精品秋霞免费鲁丝片| av女优亚洲男人天堂| 亚洲国产日韩一区二区| 亚洲国产av新网站| 亚洲一区中文字幕在线| 国产午夜精品一二区理论片| 亚洲精品视频女| 国产精品免费大片| 一级黄片播放器| 欧美+日韩+精品| 精品一区二区三区四区五区乱码 | 亚洲第一av免费看| 18在线观看网站| 美女福利国产在线| 久久久精品免费免费高清| 久久久精品国产亚洲av高清涩受| 国产精品久久久久久精品电影小说| 欧美+日韩+精品| 各种免费的搞黄视频| 精品国产超薄肉色丝袜足j| 中文字幕人妻丝袜制服| 日本爱情动作片www.在线观看| 在线观看人妻少妇| 欧美精品一区二区免费开放| 亚洲国产最新在线播放| 尾随美女入室| 中文精品一卡2卡3卡4更新| 欧美日韩亚洲高清精品| 中文字幕人妻丝袜制服| 欧美另类一区| 纯流量卡能插随身wifi吗| 制服丝袜香蕉在线| 人人澡人人妻人| 午夜福利视频在线观看免费| 香蕉国产在线看| 2021少妇久久久久久久久久久| 日韩在线高清观看一区二区三区| 精品久久久久久电影网| 麻豆乱淫一区二区| 一级毛片黄色毛片免费观看视频| 在线观看美女被高潮喷水网站| 在线观看一区二区三区激情| 亚洲色图 男人天堂 中文字幕| 国产在视频线精品| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 欧美中文综合在线视频| 国产深夜福利视频在线观看| 日产精品乱码卡一卡2卡三| 一本久久精品| 亚洲精品美女久久久久99蜜臀 | 国产免费又黄又爽又色| 午夜av观看不卡| 国产精品久久久久久av不卡| 伊人久久国产一区二区| 精品人妻偷拍中文字幕| 美女xxoo啪啪120秒动态图| 女人被躁到高潮嗷嗷叫费观| 色网站视频免费| 男女高潮啪啪啪动态图| 国产成人欧美| 中文字幕人妻丝袜制服| 国产视频首页在线观看| 如何舔出高潮| 叶爱在线成人免费视频播放| 国产精品三级大全| 国产精品久久久久久精品古装| 久久女婷五月综合色啪小说| 国产男人的电影天堂91| 少妇人妻精品综合一区二区| 熟女电影av网| 高清av免费在线| 麻豆乱淫一区二区| 美女高潮到喷水免费观看| 日韩大片免费观看网站| 一二三四在线观看免费中文在| 十八禁网站网址无遮挡| 免费不卡的大黄色大毛片视频在线观看| 国产黄色视频一区二区在线观看| 91精品三级在线观看| 日本vs欧美在线观看视频| 国产片特级美女逼逼视频| 亚洲男人天堂网一区| 精品一区二区免费观看| 最近2019中文字幕mv第一页| 国产精品免费视频内射| 中文字幕精品免费在线观看视频| 美女视频免费永久观看网站| 看免费成人av毛片| 久久久久国产一级毛片高清牌| 日韩一区二区三区影片| 久久精品人人爽人人爽视色| 国产精品免费大片| 如何舔出高潮| 久久精品国产亚洲av高清一级| 黑人欧美特级aaaaaa片| 极品少妇高潮喷水抽搐| 最新的欧美精品一区二区| 毛片一级片免费看久久久久| 美女中出高潮动态图| 午夜av观看不卡| 亚洲精品国产av成人精品| 蜜桃在线观看..| 久久久精品国产亚洲av高清涩受| av女优亚洲男人天堂| 亚洲精品美女久久av网站| 亚洲精品乱久久久久久| 人妻系列 视频| 色网站视频免费| 欧美另类一区| 久久精品国产亚洲av天美| 美女大奶头黄色视频| 91在线精品国自产拍蜜月| 一个人免费看片子| 亚洲精品久久久久久婷婷小说| 99久久综合免费| 少妇的逼水好多| 视频区图区小说| 秋霞在线观看毛片| 黑丝袜美女国产一区| 搡老乐熟女国产| 亚洲国产毛片av蜜桃av| 一区在线观看完整版| 久久久久精品性色| 亚洲成人av在线免费| 亚洲av电影在线观看一区二区三区| 亚洲精品一区蜜桃| 亚洲av国产av综合av卡| 一级毛片我不卡| 亚洲精品美女久久久久99蜜臀 | 晚上一个人看的免费电影| 寂寞人妻少妇视频99o| 亚洲精品久久午夜乱码| 精品国产一区二区三区四区第35| 欧美人与性动交α欧美精品济南到 | 性少妇av在线| 黄网站色视频无遮挡免费观看| 国产精品av久久久久免费| 国产一区二区 视频在线| 午夜福利视频在线观看免费| 91aial.com中文字幕在线观看| 久久国内精品自在自线图片| 国产精品二区激情视频| 一级片'在线观看视频| 久久久久久人妻| 午夜日本视频在线| 亚洲精品国产av成人精品| 精品人妻熟女毛片av久久网站| 高清av免费在线| 亚洲av在线观看美女高潮| 超色免费av| 国语对白做爰xxxⅹ性视频网站| 久久av网站| 一边亲一边摸免费视频| 韩国精品一区二区三区| 在线观看免费日韩欧美大片| 久久这里只有精品19| 国产精品女同一区二区软件| 精品视频人人做人人爽| 菩萨蛮人人尽说江南好唐韦庄| 三级国产精品片| 伊人亚洲综合成人网| 日韩欧美精品免费久久| 欧美精品高潮呻吟av久久| av福利片在线| 国产精品一二三区在线看| 亚洲国产精品一区二区三区在线| 免费观看无遮挡的男女| 欧美老熟妇乱子伦牲交| 一级片'在线观看视频| 欧美av亚洲av综合av国产av | 欧美精品人与动牲交sv欧美| 亚洲欧美一区二区三区久久| 国产探花极品一区二区| 一级毛片电影观看| 欧美在线黄色| 最近最新中文字幕免费大全7| 欧美成人精品欧美一级黄| 日韩av在线免费看完整版不卡| 国语对白做爰xxxⅹ性视频网站| 亚洲在久久综合| 伦理电影免费视频| 日韩中文字幕视频在线看片| 精品99又大又爽又粗少妇毛片| 国产精品偷伦视频观看了| 97人妻天天添夜夜摸| 永久网站在线| 一区二区三区乱码不卡18| 日韩熟女老妇一区二区性免费视频| 成人18禁高潮啪啪吃奶动态图| 日韩电影二区| xxx大片免费视频| 精品人妻一区二区三区麻豆| 日韩伦理黄色片| 丝袜在线中文字幕| 青草久久国产| 少妇熟女欧美另类| 日韩av不卡免费在线播放| 国产精品嫩草影院av在线观看| 欧美日韩综合久久久久久| 丝袜美足系列| 不卡视频在线观看欧美| 成人漫画全彩无遮挡| 欧美老熟妇乱子伦牲交| 男女国产视频网站| h视频一区二区三区| 精品午夜福利在线看| 少妇精品久久久久久久| 女人精品久久久久毛片| 亚洲伊人色综图| 国产精品秋霞免费鲁丝片| 女的被弄到高潮叫床怎么办| 波多野结衣一区麻豆| 亚洲精品中文字幕在线视频| av在线观看视频网站免费| 国产乱来视频区| 欧美亚洲 丝袜 人妻 在线| 一级黄片播放器| 日韩一本色道免费dvd| 久久久久精品性色| 国产精品 国内视频| 水蜜桃什么品种好| 一级,二级,三级黄色视频| 人妻少妇偷人精品九色| 成人毛片a级毛片在线播放| 国产 精品1| 亚洲av国产av综合av卡| 建设人人有责人人尽责人人享有的| 国产极品天堂在线| 黄色 视频免费看| 日日撸夜夜添| 亚洲国产欧美日韩在线播放| 国产成人免费观看mmmm| 两个人免费观看高清视频| 青春草亚洲视频在线观看| 国产精品一区二区在线观看99| 国产精品久久久久成人av| 日韩av不卡免费在线播放| 蜜桃国产av成人99| 超色免费av| 久久精品国产鲁丝片午夜精品| 熟女少妇亚洲综合色aaa.| 亚洲国产毛片av蜜桃av| 大码成人一级视频| 欧美在线黄色| 黄色怎么调成土黄色| 亚洲人成77777在线视频| 爱豆传媒免费全集在线观看| 日本91视频免费播放| 2018国产大陆天天弄谢| 亚洲第一青青草原| 亚洲国产欧美网| 电影成人av| 女的被弄到高潮叫床怎么办| 日韩三级伦理在线观看| 亚洲第一青青草原| 日本av手机在线免费观看| 久久鲁丝午夜福利片| 丰满少妇做爰视频| 日韩av在线免费看完整版不卡|