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

    Coupled-generalized nonlinear Schr¨odinger equations solved by adaptive step-size methods in interaction picture

    2023-03-13 09:19:14LeiChen陳磊PanLi李磐HeShanLiu劉河山JinYu余錦ChangJunKe柯常軍andZiRenLuo羅子人
    Chinese Physics B 2023年2期
    關(guān)鍵詞:陳磊河山

    Lei Chen(陳磊) Pan Li(李磐) He-Shan Liu(劉河山) Jin Yu(余錦)Chang-Jun Ke(柯常軍) and Zi-Ren Luo(羅子人)

    1Aerospace Information Research Institute,Chinese Academy of Sciences,Beijing 100094,China

    2National Microgravity Laboratory,Institute of Mechanics,Chinese Academy of Sciences,Beijing 100190,China

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

    Keywords: nonlinear optics, optical propagation in nonlinear media, coupled-generalized nonlinear Schr¨odinger equations(C-GNLSE),adaptive step-size methods

    1.Introduction

    Coupling an ultrashort laser pulses into an optical fiber,a wealth of nonlinear effects will take place due to the dispersion and the nonlinear effect of the fiber.The phenomenon has been widely used in optical fiber communication, ultrafast optical,supercontinuum generation,optical coherence tomography,etc.[1-4]The propagation of the low power ultrashort pulses in the fiber could be described by a mathematical model of the nonlinear Schr¨odinger equation (NLSE) which contains the group velocity dispersion (GVD) and self-phase modulation (SPM) terms.[2]To evaluate the high peak power femtosecond pulses, the generalized nonlinear Schr¨odinger equation (GNLSE) with the high order dispersion and nonlinear terms was adopted.[5]Although the NLSE could be analytically solved by the inverse scattering and self-similar method,[6,7]the GNLSE can only be calculated numerically.A common numerical solution to the GNLSE has been obtained by the split-step Fourier method (SSFM).[2]In the SSFM scheme, the dispersion and nonlinearity were integrated respectively in each step and the accuracy of the global error was of the second-order.[8]While the accuracy could be improved by the symmetric split-step Fourier method(S-SSFM)[9-13]or high order split-step schemes such as the fourth-order scheme of Blow and Wood,[14]the global accuracy was not better than that of the nonlinear step integration.In order to further improve the accuracy,a number of methods have been developed:A Runge-Kutta in the interaction picture(RK4IP)method was extended to solve the GNLSE and stimulate the pulse propagation and supercontinuum generation in the optic fiber by Hult.[15]The RK4IP exhibits a fifth order local accuracy and high calculation efficiency in the fixed step methods.The conservation quantity error adaptive step-control algorithm based on RK4IP(RK4IP-CQE)was introduced to solve the GNLSE by Heidt in 2009.[16]The RK4IP-CQE,an effective and accurate numerical method of GNLSE solution,reduced more than 50% computational time than the local error adaptive stepcontrol method.[17]Furtherly, the RK4IP-CQE in frequency domain was also introduced to solve the GNLSE by Riezniket al.[18]and the numerical results seems impressive.

    Recently multimode optical fibers and micron waveguides have reemerged as a viable platform to observe the novel linear and nonlinear physical phenomena.[19-23]The additional spatial degrees of freedom of the new fibers and waveguides offer further opportunities to investigate the interesting phenomena and processes such as the vector soliton fission,the vector modulation instability,the intermodal modulational instability, and the soliton capture appeared.[24-33]Understanding the ultrashort pulse propagation dynamics mechanism behind these complex physics phenomena and processes needs to solve the numerical modeling of the so-called coupled-generalized nonlinear Schrodinger equations (C-GNLSE) or the multimode generalized nonlinear Schr¨odinger equations(MM-GNLSE).

    Although the MM-GNLSE can be easily handled with the SSFM and S-SSFM, the accuracies andyefficiencies of the two methods need to be enhanced.To obtain accurate and effective numerical simulations,the fixed step method of RK4IP is extend to solve MM-GNLSE.[34-36]In order to improve the calculation speed of FFT, the GPU-accelerated method has been adopted.[37]Practically,an algorithm with various adaptive step sizes such as RK4IP-CQE or local error adaptive stepcontrol method used in GNLSE simulation may be more effective and useful than the fixed step method.However, to our knowledge, the adaptive step-size algorithm has not been extended to the solution of MM-GNSLE.

    In this work we extend the adaptive step-size algorithm to C-GNLSE to verify the effectiveness.By mapping the CGNLSE in the normal picture into the interaction picture in frequency domain and converting the coupled equations into a vector equation, the RK4IP-CQE and the local error adaptive step-control method are introduced to solve C-GNLSE.The two adaptive step-control algorithms are used to solve CGNLSE to simulate the SC generation in high birefringence PCF.The simulation results are the same as those calculated by fixed step algorithm RK4IP, which proves the accuracies of the two algorithms.The calculation efficiency of RK4IPLEM and RK4IP-CQE are displayed and the large difference in computational time between the two algorithms at the same global error is explained.

    2.Coupled generalized nonlinear Schro¨dinger equations

    To simplify our numerical simulations,RK4IP-CQE and RK4IP-LEM are used for solving the two-dimensional CGNLSE in instead of MM-GNLSE.The process and the methods can be fully extended to solve then-dimensional MMGNLSE withndimensions.The typical coupled generalized nonlinear Schr¨odinger equations (C-GNLSE) can be expressed as[38-40]

    with

    whereAnis the pulse envelope for the polarizationn, the asterisk denotes the complex conjugate, the timeT=tz(β11+β12)/2 is in a reference frame moving at a group velocity,βmn=?βn/?ω|ω=ωnis them-th term of the Taylor series expansion for the propagation constantβn(ω),δβ1=(β11-β12)/2, Δβ=β01-β02,γn=n2ω0/cAeff,n,cis the speed of light in vacuum,Aeff,nis the effective area for the polarizationn,andω0is the carrier frequency.The express of shock term isτSHOCK,n=1/ω0.R1(T),R2(T),andR3(T)are the response functions of the fiber which is expressed as

    wherefR=0.18 is the Raman response contribution to the Kerr effect.f1(t) andf3(t) are related to the parallel and orthogonal Raman gain respectively, and can be measured experimentally.[41,42]In the simulations, the expressions off1(t),f2(t), andf3(t) aref1= (τ21+τ22)/τ1τ22×exp(-t/τ2)sin(t/τ1),f2=f1- 2f3, andf3= [(2τb-t)/τ2b)]exp(-t/τb) respectively and the values ofτ1,τ2, andτbare 12.2 fs, 32 fs, and 96 fs respectively.The C-GNLSE shown in Eq.(1) can be converted to the form indicated in Refs.[31,32,43] if the vertical and the parallel Raman gain functions are omitted.

    Equation(1)can be translated to the frequency domain by Fourier transform as follows:

    3.Algorithm

    3.1.Fourth-order Runge-Kutta in interaction picture method

    Equation(5)can be changed into the vector form,

    The final solution of C-GNSLE can be obtained by the integration of each step using the iterative scheme of RK4IP at the fifth-order local accuracy.However, the error of RK4IP method in integration step cannot be predetermined unless the step is small enough.In order to control the error and change the step within the error in the stimulations, the local error adaptive step method and conservation quantity error adaptive step method based on RK4IP will be introduced.

    3.2.RK4IP-LEM

    For the C-GNSLE in Eq.(8), if the complex field is discretized into the frequency grid points and an integration method with RK4IP in Eq.(9)is uesed,there exists a constant for each grid point so that the calculated field can be expressed as

    where ?Ature(z+h,ω) is an exact solution.Here,η=5 for RK4IP scheme,for it exhibits a five-order local error.The relative local amplitude error is now defined as

    3.3.RK4IP-CQE

    The optical photon numberPduring the propagation can be given by

    where

    4.Results and discussion

    In this section, the performance of adaptive step algorithms described in the above section is compared and discussed.

    To prove the accuracy of the methods for C-GNLSE,a typical example of the supercontinuum generation in high birefringence photonic crystal fiber(PCF)is first simulated by using a constant step size of RK4IP.The input pulse is a hyperbolic secant with a full width at half maximum (FWHM)durationTFWHM=50 fs.The peak power of the input pulse is 20 kW, and the center wavelength is 680 nm.The angle between the polarization and the fast axis of the high birefringence PCF isπ/4.The length of high birefringence PCF is 0.1 m,and the Taylor expansion coefficients for the dispersion curve are taken from Martinset al.[39]The nonlinear coefficient isγ1=γ2=0.045 W-1·m-1.In the stimulation,the time window is 5 ps and discretized into 213grids.

    Figure 1 illustrates a temporal and spectral evolution of the supercontinuum generation process over 0.1-m length of high birefringence PCF with the step size of 40 μm in the stimulation using the constant step-size method of RK4IP.A logarithmic density scale is used which is truncated at-80 dB relative to the maximum value.As shown in Figs.1(a)-1(d),the orthogonally polarized pulses traveling at different group velocities in the slow axial direction and the fast axial direction are completely separated in time after a few-mm propagation distance.After the orthogonally polarized pulses transmit apart,both of them break into a series of pulses at about 1.5 cm and 2 cm far in the slow axial direction and fast axial direction respectively, which is known as the vector soliton fission.[30]The fundamental solitons emerge one by one in the slow axial direction and the fast axial direction and subsequently shift to longer wavelengths due to the intrapulse Raman scattering.Therefore, the energy is transferred to a narrow band resonance in the normal GVD regime, associated with the emergence of a dispersive wave.The spectra of the different red moving solitons and dispersion waves are formed as an octave supercontinuum.The simulations of the supercontinuum generated in high birefringence PCF are similar to the results obtain by Martinset al.[39]

    Fig.1.Numerical simulation of supercontinuum generation in 0.1-mhigh birefringence PCF.(a) Temporal evolution and (b) spectral evolution of pulse along slow axis; (c) temporal evolution and (d) spectral evolution of pulse along fast axis, with retarded time frame of the reference travelling at envelope group velocity of input pulse used in panels(a)and(b).

    To exhibit the errors changing with propagation distance,fgiure 2 shows the approximate local error and the relative photon number error between the values of the two consecutive computational steps with the pulse transmitting in the high birefringence PCF.The calculation is run with RK4IP method in the vector form in steps of 40 μm.The approximated local error and the relative photon number error are calculated with Eqs.(11)and(16).As shown in Figs.2(a)and 2(b),the approximated local error curve is similar to the relative photon number error curve, where they both have a large error peak located in a range from 0 m to 0.01 m and a small error peak located at about 0.02 m.To further understand the details of the error variance, theyaxes of the approximated local error curve and relative photon number error curve are enlarged as indicated in the insets, respectively, in Figs.2(a)and 2(b).The details of two curves are also similar to each other and they both fluctuate heavily during the propagation in a distance of 0 m-0.05 m; while the propagating distance is more than 0.05 m,the two error curves rise monotonically.

    The oscillations of the two error curves ranging from 0 m to 0.05 m are mainly due to the vector soliton fission process.Many different frequency photons are created and annihilated in the vector soliton fission process,which makes the spectrum of the pulse change and the photon number error increase.At the same time, the time waveform of the pulse breaks into a series of solitons, which increases the local error.When the propagation distance is larger than 0.05 m, the vector soliton fission and photons creation and annihilation become weak.Therefore the photon number error and local error vary slowly.

    Fig.2.Plots of propagation distance-dependent error estimations between two consecutive computational steps for (a) approximate local error(Eq.(11)),and(b)relative photon number error(Eq.(16))in a constant step of 40 μm, with inset showing enlarged section to compare the small scale errors,and error estimations depicted with the same x axis to facilitate graph comparison.

    From the above analysis, it can be concluded that when the waveforms in the time domain and in frequency domain change dramatically,the approximate local error curve and relative photon number error become large,so,the step should be reduced to produce a small error;otherwise,when the approximated local error and relative photon number error become small,step should be increased to reduce the simulation time.

    Supercontinuum generation in a high-birefringence photonic crystal fiber (PCF) is also simulated with the same parameters by RK4IP-LEM with a relative local errorηG=10-6and RK4IP-CQE with a relative photon number errorηP=10-14.Figure 3 shows the calculation results.The output spectra calculated by three different methods are completely the same,which means that RK4IP-LEM and RK4IP-CQE are reliable.

    Fig.3.Slow and fast axis output spectra of high birefringence PCF calculation by RK4IP,RK4IP-LEM,and RK4IP-CQE methods.

    To make a comparison of accuracy and efficiency between the RK4IP-LEM method and the RK4IP-CQE method,the computational time value with the global error of the twoadaptive step-control algorithm is shown in Fig.4.The global average error in Fig.4 is defined as

    where the complex field ?Acalis calculated by the RK4IP-LEM method and the RK4IP-CQE method, ?Aaccis calculated by the RK4IP in steps of 0.1μm.

    It can be seen from Fig.4 that the computational time taken by each of the RK4IP-LEM method and the RK4IPCQE method generally decreases with the global average error increasing.However,the computational time taken by the RK4IP-LEM method is much more than that by the RK4IPCQE method,which is 30 times more than that by the RK4IPCQE method when the global average error is less than 10-7,and is about 20 times more than that by the RK4IP-CQE method when the global average error is more than 10-7.

    Fig.4.Computational time taken by RK4IP-CQE and RK4IP-LEM against global average error, normalized by the time required to evaluate 103 FFT, for the supercontinuum generation process in highbirefringence PCF.

    The RK4IP-LEM method uses one more FFTs than the RK4IP-CQE at each step as indicated in the above section.If the step sizes of two methods in each iteration process are the same, the computational time of RK4IP-LEM will be twice longer than that of the RK4IP-CQE method.In order to explain the difference in computational time between the RK4IPLEM method and the RK4IP-CQE method, the error estimations between the two consecutive steps for the approximate local error and relative photon number error are calculated by the RK4IP method in different step sizes,furthermore,the valuations of the step size and the error between two consecutive steps in calculation process of the RK4IP-LEM method and the RK4IP-CQE method under different error limits are also calculated.The results are shown in Figs.5-7,respectively.

    Figure 5 shows the approximate local error and relative photon number error between two consecutive steps calculated by the RK4IP method in steps of 10μm,20μm,40μm,and 80 μm.As shown in Fig.5(a), the approximate local error curves are similar when the step size is 40 μm and 80 μm,respectively.They both fluctuate heavily in a propagation distance range of 0 m-0.05 m.When the pulse propagating distance is more than 0.05 m,the two error curves rise monotonically.The reason of the fluctuations can be found in Fig.2.While the step size decreases to 10μm and 20μm,the fluctuation between 0 m-0.05 m disappear,which is different from that in steps of 40 μm and 80 μm.As shown in Fig.5(b),except for different error amplitudes,the relative photon number curves are all nearly the same in steps of 80 μm, 40 μm,20 μm, and 10 μm.They all fluctuate heavily in the propagation distance of 0 m-0.05 m and rise monotonically for the propagation distance large than 0.05 m.

    Fig.5.Variations of (a) approximated local error (Eq.(11)) and (b)relative photon number error(Eq.(17))with different propagation distances,calculated by RK4IP method in steps of 10μm,20μm,40μm,80μm,with y axis being logarithmic.

    Figure 6 shows the variations of step size with the propagation distance in the process by using the RK4IP-LEM method and the RK4IP-CQE method under different approximate local error limitsηGand the relative photon number error limitsηPrespectively.In Fig.6(a), except small oscillations near 0.01 m and 0.02 m,the step size increases until the propagation distance reaches 0.05 m at the approximate local error limitηG=10-5.The step size decreases monotonically after 0.05 m.The small oscillations near 0.01 m and 0.02 m disappear gradually and the step size increases monotonically before 0.05 m and then decreases monotonically as theηGdecreases to 10-6,10-7,and 10-8.The smaller theηG,the more gently the change of the step size is.It is shown in Fig.6(b)that the step size curves with different relative photon number limitsηP=10-12; 10-13; 10-14, 10-15all oscillate heavily in a propagation distance range of 0 m-0.03 m; the step size increases monotonically in 0.03 m-0.05 m and the step size decreases monotonically in 0.05 m-0.1 m.

    Fig.6.Step sizes versus propagation distances obtained by(a)RK4IPLEM)and(b)RK4IP-CQE under different error limits.

    Fig.7.(a) Approximated local error of RK4IP-LEM and (b) relative photon number error of RK4IP-CQE between consecutive computational steps under different error limits.

    Figure 7 is the approximated local error(Fig.7(a))and the relative photon number error(Fig.7(b))between the consecutive computational steps varying with the propagation distance in solving the C-GNLSE by the RK4IP-LEM method and the RK4IP-CQE method with different error limits.In Fig.7(a)the approximate local errors are all under the presetting error limits atηGequating 10-5,10-6,10-7,and 10-8.The fluctuations of the approximate local error curves become gentle as the change time of the step (ηG) decreases.In Fig.7(b), except within the propagation distance between 0.01 m-0.02 m,the relative photon number errors are all under the presetting error limits atηPequating 10-12, 10-13, 10-14, and 10-15.The fluctuations of the relative photon number error are nearly the same as those as theηPdecreases, which means that the change of the step is independent of the value ofηP.

    From the above analyses, the difference between computational time taken by the RK4IP-LEM method and the computational time taken by the RK4IP-CQE method at the same global average error level in Fig.4 can be qualitatively explained by the results of Figs.5-7.Owing to the nonconserved qualitatively approximated local error,the approximated local error curves under different step sizes are not similar(see Fig.5(a)), which means that convergence rate of the approximated local error is much different from the error limitηGat every step of the propagation distance.The different convergence rates make the change of step size and approximated local error in the process of RK4IP-LEM different(see Figs.6(a) and 7(a)).The non-similarity of the change of the step size makes global average error of RK4IP-LEM not uniformly converge.However, the relative photon number error curves under different step sizes are similar because the relative photon number error is conserved quantity(see Fig.5(b)),and the change of step size and the relative photon number in the process of RK4IP-CQE method under different error limits are similar (see Figs.6(b) and 7(b)).The similarity of the changing of the step size makes the global average error of RK4IP-CQE method converge uniformly.The difference between the convergences of the relative photon number error and the approximated local error induces the large difference in computational time at the same global average error in Fig.4.

    5.Conclusions

    By mapping the C-GNLSE in the normal picture into the interaction picture in the frequency domain, the conservation quantity error adaptive step-control method and the local error adaptive step-control method are developed based on a vector form of the fourth-order Runge-Kutta in interaction picture.To prove the efficiency of the adaptive step-control methods, the two adaptive step-control methods and the RK4IP method are used to simulate the SC generation in the high birefringence PCF.The calculation accuracy and efficiency for each of these two adaptive step-control methods are discussed.At the same global average error, the computational time of RK4IP-CQE has been improved 20 times compared with that of RK4IP-LEM due to the convergences of the relative photon number error and the approximated local error.The methods will be useful for simulating the vector pulses transmission and the supercontinuum generation in the nonlinear fiber and waveguides, the pulse propagating in the multimode optical fiber,and the interaction between different pulses.

    Appendix A

    LetAnbe the pulse envelope for the polarizationn, then the photon numberPduring propagation will be defined as

    Calculating each integral term in Eq.(A3),the following equation are obtained:

    LetΩ=ω1,Ω1=ω,Ω2=ω2+ω1-ωand use the expression ofRn(Ω1-Ω)=R?n(Ω-Ω1), then equations(A4)and(A5)will change into Eqs.(A8)and(A9),i.e.,

    Substituting Eqs.(A6)-(A9)into Eq.(A3),then the expression of?P/?z=0 is obtained.

    Acknowledgement

    Project supported by the National Key Research and Development Program of China (Grant Nos.2021YFC2201803 and 2020YFC2200104).

    猜你喜歡
    陳磊河山
    容易到摳腳的相對論
    容易到摳腳的相對論
    陳磊
    河山
    中國寶玉石(2021年5期)2021-11-18 07:34:50
    楚漢之爭(七)
    楚 漢 之 爭(八)
    故鄉(xiāng)
    北方音樂(2020年7期)2020-06-01 07:26:13
    情系河山
    中華詩詞(2019年12期)2019-09-21 08:53:06
    情系河山
    中華詩詞(2018年5期)2018-11-22 06:46:10
    情系河山
    中華詩詞(2018年1期)2018-06-26 08:46:30
    日本 欧美在线| 午夜激情av网站| 日韩有码中文字幕| 97人妻精品一区二区三区麻豆 | 国产精品久久电影中文字幕| 后天国语完整版免费观看| 国产午夜精品久久久久久| 久久99热这里只有精品18| 韩国精品一区二区三区| 成人国产综合亚洲| 久久精品国产综合久久久| 精品一区二区三区av网在线观看| 国产单亲对白刺激| 黄频高清免费视频| 亚洲av五月六月丁香网| 中出人妻视频一区二区| 成人免费观看视频高清| 午夜影院日韩av| 超碰成人久久| 国内久久婷婷六月综合欲色啪| 久9热在线精品视频| 琪琪午夜伦伦电影理论片6080| 久久久国产欧美日韩av| 老熟妇乱子伦视频在线观看| 怎么达到女性高潮| 精品午夜福利视频在线观看一区| 久久午夜亚洲精品久久| 男人的好看免费观看在线视频 | 中文字幕高清在线视频| 老司机午夜福利在线观看视频| 香蕉国产在线看| 桃红色精品国产亚洲av| 天堂√8在线中文| 日本免费一区二区三区高清不卡| 757午夜福利合集在线观看| 成人欧美大片| 午夜免费鲁丝| 午夜福利成人在线免费观看| 精品无人区乱码1区二区| 制服诱惑二区| 一区二区三区精品91| 亚洲黑人精品在线| 色播亚洲综合网| 免费看美女性在线毛片视频| 亚洲熟妇熟女久久| 国产午夜精品久久久久久| 亚洲最大成人中文| 99热这里只有精品一区 | 精品无人区乱码1区二区| 亚洲aⅴ乱码一区二区在线播放 | 精品人妻1区二区| 黑人巨大精品欧美一区二区mp4| av电影中文网址| 亚洲av片天天在线观看| 中文字幕久久专区| 国产aⅴ精品一区二区三区波| 午夜久久久久精精品| 精品国产乱子伦一区二区三区| 精品久久久久久久久久久久久 | 久久 成人 亚洲| 国产成人精品无人区| 他把我摸到了高潮在线观看| 一区二区日韩欧美中文字幕| 成人av一区二区三区在线看| 亚洲一区二区三区不卡视频| 日本一本二区三区精品| e午夜精品久久久久久久| 两性夫妻黄色片| 亚洲精品中文字幕一二三四区| 亚洲熟妇熟女久久| 精品国产美女av久久久久小说| 日韩 欧美 亚洲 中文字幕| 国产精品久久久久久亚洲av鲁大| 精品国产亚洲在线| 成人特级黄色片久久久久久久| 国产熟女xx| 国产精品野战在线观看| 在线播放国产精品三级| 两性夫妻黄色片| 后天国语完整版免费观看| 99在线视频只有这里精品首页| 日韩高清综合在线| 国产欧美日韩一区二区精品| 亚洲欧美一区二区三区黑人| 免费看美女性在线毛片视频| 琪琪午夜伦伦电影理论片6080| 男女床上黄色一级片免费看| 深夜精品福利| 日韩欧美一区视频在线观看| 国产一区在线观看成人免费| 精品国产一区二区三区四区第35| 巨乳人妻的诱惑在线观看| 欧美av亚洲av综合av国产av| 欧美又色又爽又黄视频| 国产亚洲精品第一综合不卡| 久久人人精品亚洲av| 久久久国产欧美日韩av| av在线播放免费不卡| 老司机在亚洲福利影院| 亚洲天堂国产精品一区在线| 亚洲第一欧美日韩一区二区三区| 俺也久久电影网| 啦啦啦 在线观看视频| 欧美+亚洲+日韩+国产| 午夜激情av网站| 久久国产精品影院| 欧美绝顶高潮抽搐喷水| 老司机深夜福利视频在线观看| 少妇熟女aⅴ在线视频| 亚洲av五月六月丁香网| 婷婷精品国产亚洲av| 国产精品98久久久久久宅男小说| 日本黄色视频三级网站网址| 日本熟妇午夜| 免费观看精品视频网站| 听说在线观看完整版免费高清| 色在线成人网| 日本免费一区二区三区高清不卡| 国产精品亚洲一级av第二区| 国产精品国产高清国产av| 精品卡一卡二卡四卡免费| 一级毛片精品| 淫妇啪啪啪对白视频| 精品久久久久久久久久久久久 | 我的亚洲天堂| 正在播放国产对白刺激| 一本精品99久久精品77| 制服人妻中文乱码| 免费在线观看日本一区| 亚洲人成网站高清观看| 我的亚洲天堂| 丁香六月欧美| 人人妻人人看人人澡| 日韩国内少妇激情av| 国产精品99久久99久久久不卡| 亚洲成人久久性| 亚洲欧美精品综合一区二区三区| 天天一区二区日本电影三级| 色综合婷婷激情| 老司机深夜福利视频在线观看| 曰老女人黄片| 超碰成人久久| 精品不卡国产一区二区三区| 丝袜在线中文字幕| 好男人电影高清在线观看| 亚洲美女黄片视频| 男女做爰动态图高潮gif福利片| 亚洲成人国产一区在线观看| 亚洲精品粉嫩美女一区| 欧美日韩亚洲综合一区二区三区_| 亚洲av日韩精品久久久久久密| 亚洲av熟女| 亚洲五月婷婷丁香| 亚洲无线在线观看| 日韩精品免费视频一区二区三区| 亚洲av电影在线进入| 亚洲成av片中文字幕在线观看| 精品国产一区二区三区四区第35| 欧美日韩黄片免| 白带黄色成豆腐渣| 色婷婷久久久亚洲欧美| 免费无遮挡裸体视频| 日本黄色视频三级网站网址| 伊人久久大香线蕉亚洲五| 波多野结衣av一区二区av| 久久精品aⅴ一区二区三区四区| 国产精品免费一区二区三区在线| 国产精品日韩av在线免费观看| 日韩av在线大香蕉| 国产色视频综合| 成人手机av| 亚洲色图 男人天堂 中文字幕| 99在线视频只有这里精品首页| 香蕉丝袜av| 丝袜美腿诱惑在线| 亚洲欧美一区二区三区黑人| 亚洲一码二码三码区别大吗| 日本 欧美在线| 俺也久久电影网| 99在线视频只有这里精品首页| 一进一出抽搐动态| 狂野欧美激情性xxxx| 国内久久婷婷六月综合欲色啪| 日日爽夜夜爽网站| 一本综合久久免费| 一区二区三区激情视频| 国产av在哪里看| 欧美黑人精品巨大| 美女高潮到喷水免费观看| 18美女黄网站色大片免费观看| 最近最新免费中文字幕在线| 男人操女人黄网站| 国产精品美女特级片免费视频播放器 | 中文字幕最新亚洲高清| 91老司机精品| 亚洲自偷自拍图片 自拍| 婷婷亚洲欧美| 韩国精品一区二区三区| 国产蜜桃级精品一区二区三区| 琪琪午夜伦伦电影理论片6080| 免费在线观看完整版高清| videosex国产| 亚洲av第一区精品v没综合| cao死你这个sao货| 亚洲精品久久国产高清桃花| 99riav亚洲国产免费| 久9热在线精品视频| 亚洲免费av在线视频| 成人18禁高潮啪啪吃奶动态图| 黑人操中国人逼视频| 不卡一级毛片| 老熟妇仑乱视频hdxx| 亚洲自拍偷在线| 久久中文看片网| 91av网站免费观看| 亚洲五月婷婷丁香| 特大巨黑吊av在线直播 | 婷婷精品国产亚洲av在线| 中文字幕人成人乱码亚洲影| 亚洲熟妇中文字幕五十中出| av视频在线观看入口| 欧美日韩一级在线毛片| 国产精品久久电影中文字幕| 人人妻人人看人人澡| 亚洲人成网站在线播放欧美日韩| 欧美日本亚洲视频在线播放| 琪琪午夜伦伦电影理论片6080| 午夜福利一区二区在线看| 成人亚洲精品一区在线观看| 国产av不卡久久| av有码第一页| 日韩视频一区二区在线观看| 国产黄a三级三级三级人| 国产黄色小视频在线观看| 精品久久久久久久毛片微露脸| 精品久久久久久久毛片微露脸| a在线观看视频网站| 18禁黄网站禁片午夜丰满| 亚洲第一电影网av| 国产伦在线观看视频一区| 日韩中文字幕欧美一区二区| 久久国产乱子伦精品免费另类| 国产欧美日韩精品亚洲av| 啪啪无遮挡十八禁网站| 欧美日韩福利视频一区二区| 国产aⅴ精品一区二区三区波| 色播在线永久视频| 日本黄色视频三级网站网址| 波多野结衣高清无吗| 国产精品久久久av美女十八| 2021天堂中文幕一二区在线观 | 欧美乱妇无乱码| 狂野欧美激情性xxxx| 黄色a级毛片大全视频| 欧美成人午夜精品| 88av欧美| 国产精华一区二区三区| 黄色片一级片一级黄色片| 午夜成年电影在线免费观看| 精品高清国产在线一区| 国产免费男女视频| 欧美三级亚洲精品| 国产精品香港三级国产av潘金莲| 天堂√8在线中文| 桃红色精品国产亚洲av| 国产久久久一区二区三区| 无限看片的www在线观看| www.自偷自拍.com| www.999成人在线观看| 男人的好看免费观看在线视频 | 最近最新中文字幕大全免费视频| 久久这里只有精品19| 欧美一级a爱片免费观看看 | 免费观看精品视频网站| 亚洲五月色婷婷综合| 村上凉子中文字幕在线| 亚洲中文av在线| 国产欧美日韩一区二区精品| 国产成人欧美在线观看| 亚洲国产欧美网| 黄网站色视频无遮挡免费观看| 国产一区在线观看成人免费| a级毛片a级免费在线| 亚洲七黄色美女视频| 18美女黄网站色大片免费观看| 亚洲欧美激情综合另类| 午夜久久久久精精品| 欧美 亚洲 国产 日韩一| 少妇粗大呻吟视频| 国产欧美日韩一区二区三| 国产精品久久久久久人妻精品电影| 亚洲成a人片在线一区二区| 精品一区二区三区av网在线观看| 精品少妇一区二区三区视频日本电影| av在线播放免费不卡| 成人亚洲精品一区在线观看| 最近在线观看免费完整版| √禁漫天堂资源中文www| 亚洲一区中文字幕在线| 国产成人系列免费观看| 亚洲九九香蕉| 在线看三级毛片| 一级毛片高清免费大全| ponron亚洲| 老汉色av国产亚洲站长工具| 成人手机av| xxx96com| 黄网站色视频无遮挡免费观看| 一区福利在线观看| 51午夜福利影视在线观看| 99久久精品国产亚洲精品| 看黄色毛片网站| 人人妻,人人澡人人爽秒播| 久久99热这里只有精品18| 国产欧美日韩精品亚洲av| а√天堂www在线а√下载| 中文字幕人妻熟女乱码| 精品午夜福利视频在线观看一区| 色播在线永久视频| 国语自产精品视频在线第100页| 啦啦啦免费观看视频1| 1024香蕉在线观看| 欧美zozozo另类| 亚洲avbb在线观看| 亚洲自偷自拍图片 自拍| 久久国产乱子伦精品免费另类| 亚洲九九香蕉| 久久精品国产99精品国产亚洲性色| 搡老熟女国产l中国老女人| 亚洲精品一区av在线观看| 国产一区二区激情短视频| 欧美日韩福利视频一区二区| 窝窝影院91人妻| 少妇粗大呻吟视频| 久久天堂一区二区三区四区| 亚洲人成77777在线视频| 一区二区三区精品91| www.熟女人妻精品国产| 亚洲欧美精品综合一区二区三区| av天堂在线播放| 91成年电影在线观看| 精品熟女少妇八av免费久了| 欧美久久黑人一区二区| 亚洲专区国产一区二区| 久久狼人影院| 最近最新中文字幕大全免费视频| 欧美黄色淫秽网站| 国产熟女xx| 啦啦啦免费观看视频1| 满18在线观看网站| 成人亚洲精品av一区二区| 欧美激情久久久久久爽电影| 午夜福利一区二区在线看| 免费搜索国产男女视频| 中文字幕最新亚洲高清| 夜夜看夜夜爽夜夜摸| 在线永久观看黄色视频| 俺也久久电影网| 久久国产亚洲av麻豆专区| 日韩欧美国产一区二区入口| 天天添夜夜摸| 亚洲欧美精品综合久久99| 成年版毛片免费区| 日韩欧美一区二区三区在线观看| 国产色视频综合| 精品久久久久久久人妻蜜臀av| 欧美日韩黄片免| 国产午夜精品久久久久久| 在线观看日韩欧美| 午夜日韩欧美国产| 99国产极品粉嫩在线观看| 少妇裸体淫交视频免费看高清 | 国产私拍福利视频在线观看| 禁无遮挡网站| 亚洲av中文字字幕乱码综合 | 欧美日韩一级在线毛片| 欧美性长视频在线观看| 亚洲,欧美精品.| 欧美黄色淫秽网站| 黄色成人免费大全| 黄片小视频在线播放| 国产真实乱freesex| 中文字幕高清在线视频| 午夜免费激情av| 一本精品99久久精品77| 成人av一区二区三区在线看| 亚洲狠狠婷婷综合久久图片| 亚洲色图 男人天堂 中文字幕| 日日摸夜夜添夜夜添小说| 亚洲中文字幕日韩| 亚洲av成人av| 俺也久久电影网| 亚洲在线自拍视频| 精品不卡国产一区二区三区| 宅男免费午夜| 亚洲国产精品成人综合色| 欧美中文综合在线视频| 午夜激情av网站| 久久精品影院6| aaaaa片日本免费| 欧美绝顶高潮抽搐喷水| 国产精品久久久久久精品电影 | 色婷婷久久久亚洲欧美| 91麻豆精品激情在线观看国产| 欧美人与性动交α欧美精品济南到| 波多野结衣巨乳人妻| 美女大奶头视频| 丝袜人妻中文字幕| 麻豆一二三区av精品| 午夜影院日韩av| 啪啪无遮挡十八禁网站| 亚洲免费av在线视频| 又黄又粗又硬又大视频| 一a级毛片在线观看| 国产色视频综合| 国产激情偷乱视频一区二区| 波多野结衣高清无吗| 免费女性裸体啪啪无遮挡网站| 日韩欧美在线二视频| 欧美三级亚洲精品| 女性生殖器流出的白浆| 国产精品亚洲av一区麻豆| 又黄又爽又免费观看的视频| 久久天堂一区二区三区四区| 最近在线观看免费完整版| netflix在线观看网站| 精品国产一区二区三区四区第35| 又大又爽又粗| 国产精品野战在线观看| 美女免费视频网站| 午夜福利视频1000在线观看| 99国产精品99久久久久| 天天躁狠狠躁夜夜躁狠狠躁| 中国美女看黄片| 91在线观看av| 黄色女人牲交| x7x7x7水蜜桃| 亚洲自偷自拍图片 自拍| 草草在线视频免费看| 桃红色精品国产亚洲av| 99国产精品一区二区三区| 91麻豆精品激情在线观看国产| 国产成人一区二区三区免费视频网站| 亚洲国产精品合色在线| 亚洲国产欧洲综合997久久, | 久久人妻福利社区极品人妻图片| 一级片免费观看大全| 国产av一区在线观看免费| 亚洲专区中文字幕在线| 成人午夜高清在线视频 | 亚洲国产精品sss在线观看| 18禁裸乳无遮挡免费网站照片 | 欧美日韩精品网址| 国产精品久久久久久精品电影 | 欧美不卡视频在线免费观看 | 欧美av亚洲av综合av国产av| 男人操女人黄网站| 99国产精品99久久久久| 村上凉子中文字幕在线| 99国产精品一区二区蜜桃av| 国产精品久久电影中文字幕| 成人特级黄色片久久久久久久| 免费女性裸体啪啪无遮挡网站| 日本一本二区三区精品| 免费看a级黄色片| 桃红色精品国产亚洲av| 国产熟女xx| 天堂动漫精品| 少妇的丰满在线观看| 亚洲国产欧美日韩在线播放| 国产成人精品无人区| 成人精品一区二区免费| 国产欧美日韩一区二区三| 日本撒尿小便嘘嘘汇集6| 成人特级黄色片久久久久久久| 亚洲精品国产精品久久久不卡| 日本在线视频免费播放| 久久欧美精品欧美久久欧美| 高清毛片免费观看视频网站| 亚洲午夜理论影院| av在线播放免费不卡| 最近在线观看免费完整版| 激情在线观看视频在线高清| 亚洲精品美女久久av网站| 成人三级做爰电影| 99久久久亚洲精品蜜臀av| cao死你这个sao货| 看片在线看免费视频| 欧美黑人精品巨大| 午夜激情福利司机影院| 黑人巨大精品欧美一区二区mp4| 亚洲性夜色夜夜综合| 男女那种视频在线观看| 美女免费视频网站| 操出白浆在线播放| 国产亚洲精品久久久久久毛片| 丝袜美腿诱惑在线| 精品国内亚洲2022精品成人| 国产高清激情床上av| 日本撒尿小便嘘嘘汇集6| 999久久久精品免费观看国产| 最近最新中文字幕大全免费视频| 黄色毛片三级朝国网站| 91麻豆精品激情在线观看国产| 久久草成人影院| 不卡一级毛片| 夜夜爽天天搞| 又黄又粗又硬又大视频| 国产av不卡久久| 精品久久久久久久久久免费视频| 国产精品久久久久久精品电影 | 中文字幕精品免费在线观看视频| 久久香蕉国产精品| 中文在线观看免费www的网站 | 精品国内亚洲2022精品成人| 午夜免费观看网址| www.熟女人妻精品国产| 天天躁狠狠躁夜夜躁狠狠躁| 精品熟女少妇八av免费久了| 亚洲欧美精品综合一区二区三区| 欧美日韩乱码在线| 亚洲av中文字字幕乱码综合 | 性欧美人与动物交配| 久久久久久人人人人人| 男人的好看免费观看在线视频 | 精品卡一卡二卡四卡免费| 国产一卡二卡三卡精品| 性欧美人与动物交配| 亚洲精品在线观看二区| 大香蕉久久成人网| 亚洲色图av天堂| 国产亚洲精品久久久久5区| 亚洲 欧美一区二区三区| 日韩欧美三级三区| www.熟女人妻精品国产| 丁香欧美五月| 国产亚洲欧美98| 日韩精品青青久久久久久| 一本一本综合久久| 亚洲国产高清在线一区二区三 | 久久香蕉激情| 久久久久久免费高清国产稀缺| 俄罗斯特黄特色一大片| 国内精品久久久久精免费| 长腿黑丝高跟| 日本三级黄在线观看| 男女午夜视频在线观看| 99久久精品国产亚洲精品| 精品一区二区三区视频在线观看免费| 亚洲精品中文字幕一二三四区| 欧美乱码精品一区二区三区| 88av欧美| 丰满人妻熟妇乱又伦精品不卡| 无遮挡黄片免费观看| 97超级碰碰碰精品色视频在线观看| 久久久久久人人人人人| 真人一进一出gif抽搐免费| 午夜两性在线视频| www.www免费av| 在线观看66精品国产| 99久久无色码亚洲精品果冻| 少妇熟女aⅴ在线视频| 日本免费a在线| 国产成+人综合+亚洲专区| 久久久精品欧美日韩精品| 看黄色毛片网站| 视频区欧美日本亚洲| 男女那种视频在线观看| 美女高潮喷水抽搐中文字幕| 国产人伦9x9x在线观看| 88av欧美| 国内精品久久久久精免费| 天堂动漫精品| 国产成人av激情在线播放| 88av欧美| 人妻丰满熟妇av一区二区三区| 久久精品影院6| 欧美一级a爱片免费观看看 | 国产一区在线观看成人免费| 国产成人欧美在线观看| 亚洲国产毛片av蜜桃av| www国产在线视频色| 无遮挡黄片免费观看| 久久久国产成人精品二区| 午夜福利18| 老司机午夜福利在线观看视频| 波多野结衣av一区二区av| 性欧美人与动物交配| 国产成人系列免费观看| 中出人妻视频一区二区| 国产又爽黄色视频| 老司机福利观看| 国内揄拍国产精品人妻在线 | 两性午夜刺激爽爽歪歪视频在线观看 | 欧美国产精品va在线观看不卡| 亚洲精品中文字幕一二三四区| 熟妇人妻久久中文字幕3abv| 日本a在线网址| 国产乱人伦免费视频| 在线观看午夜福利视频| 亚洲人成网站在线播放欧美日韩| 午夜免费激情av| 一进一出抽搐gif免费好疼| 99久久国产精品久久久| 人人妻,人人澡人人爽秒播| 97超级碰碰碰精品色视频在线观看| 韩国精品一区二区三区| 两人在一起打扑克的视频| 免费看a级黄色片| 高清毛片免费观看视频网站| 欧美亚洲日本最大视频资源| 一级作爱视频免费观看| 男人舔奶头视频| 成年免费大片在线观看| 波多野结衣av一区二区av|