亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        電磁波在任意磁偏角等離子體中的傳播

        2017-03-09 02:54:29張景付海洋
        電波科學(xué)學(xué)報(bào) 2017年6期

        張景 付海洋

        (復(fù)旦大學(xué) 電磁波信息科學(xué)教育部重點(diǎn)實(shí)驗(yàn)室,上海 200433)

        引 言

        等離子體廣泛存在于空間電離層壞境中, 其對電磁波的反射作用,會(huì)給測控通信信號(hào)帶來較大的反射損耗,降低測控通信信號(hào)的功率,影響飛行器和地面測控站之間測控通信信號(hào)的傳輸.對于等離子體隱身技術(shù)而言,則正是利用等離子體吸收能力強(qiáng),吸收帶寬大的特點(diǎn),完成對雷達(dá)波的吸收,使電磁波無法反射回到雷達(dá)接收機(jī). 近年來,波和等離子體的非線性相互作用在電離層加熱和黑障問題等方面逐漸得到重視[1-2]. 因此,研究電磁波在等離子體中傳播,對空間通信和空間環(huán)境感知具有重要意義.

        對于電磁波在等離子體中傳播問題的研究, Laroussi等人采用溫采爾-克勞麥斯-勃立魯英 (Wentzel-Kramers-Brillouin, WKB) 方法研究了磁化等離子體各參數(shù)對電磁波的反射、透射和吸收特性[3];Hu 等人采用散射矩陣(Scatter Matrix Method, SMM)方法分析計(jì)算了電子密度指數(shù)和拋物線分布的非均勻磁化等離子體層的反射、透射和吸收[4]; 柳文和Croft等人采用射線追蹤求解電磁波穿過電離層的傳播[5-6];文獻(xiàn)[7-13]提出了各種求解等離子體的時(shí)域有限差分(Finite Difference Time Domain, FDTD)方法等等. WKB和 SMM方法都只能求解平行分層介質(zhì)的傳播,WKB是對解析解的近似,對介質(zhì)的分布條件要求極高,必須是緩變介質(zhì),而射線追蹤法是一種基于幾何光學(xué)的高頻近似,求解頻率受限. 這幾種方法對復(fù)雜問題的求解都具有一定局限性,無法處理波和粒子的非線性現(xiàn)象.

        FDTD方法則可以求解任意等離子體分布,特別可求解波和粒子相互作用等共振現(xiàn)象.已有的FDTD方法大都只求解了電磁波傳播方向與背景磁場平行的情況,然而實(shí)際對于像電離層這種任意磁偏角的環(huán)境是比較常見的. 近期,不同學(xué)者采用梯形遞歸卷積時(shí)域有限差分(Trapezoidal Recursive Convolution, TRC-FDTD)[9],移位算子時(shí)域有限差分(Shift Operator, SO-FDTD)[12],分段線性電流密度卷積時(shí)域有限差分(Piecewise Linear Recursive Convolution, PLRC-FDTD)[14]計(jì)算了任意磁偏角的傳播,但是在計(jì)算時(shí)引入中間變量較多,占用內(nèi)存較大,計(jì)算效率低. 為簡化計(jì)算,中間變量忽略了虛部計(jì)算,且只對幅頻特性進(jìn)行了研究,并未研究相頻特性.

        本文選擇計(jì)算精度高、效率高、占用內(nèi)存少的JEC-FDTD方法[13],推導(dǎo)獲得了任意磁偏角的JEC-FDTD迭代公式. 首先計(jì)算了磁偏角θ=π/2不同頻率電磁波在等離子體傳播的反射/透射系數(shù)幅度和相位,分析了法拉第旋轉(zhuǎn)效應(yīng),驗(yàn)證算法的準(zhǔn)確性和精度.然后,分析了任意磁偏角電磁波的傳播.最后,進(jìn)一步研究了等離子體在電子回旋頻率和高頻混雜頻率處的共振吸收特性,并結(jié)合電離層加熱參數(shù),定量研究了電離層加熱高頻混雜共振吸收特性,對今后電離層加熱具有指導(dǎo)意義.

        1 任意磁偏角的JEC-FDTD算法

        對于磁化碰撞冷等離子體,假設(shè)電磁波沿+z方向傳播,外部磁場B位于yoz平面內(nèi)且與+z軸夾角為θ,忽略離子運(yùn)動(dòng),則Maxwell方程及本構(gòu)關(guān)系可寫成如下形式:

        (1)

        (2)

        (3)

        (4)

        為了便于采用FDTD迭代,令電流密度分量空間節(jié)點(diǎn)和電場分量空間節(jié)點(diǎn)位于同一節(jié)點(diǎn),電流密度分量位于半時(shí)間步,采用中心差分近似代替微分[13-14].式(4)的電流密度各分量的FDTD迭代方程如下所示:

        (5)

        (6)

        (7)

        將式(5)~(7)代入式(1)和(2)可得到電場和磁場各分量的迭代方程.

        2 數(shù)值計(jì)算和分析

        2.1 算例驗(yàn)證

        2.1.1 磁偏角θ=π/2的傳播特性

        為驗(yàn)證任意磁偏角JEC-FDTD算法的有效性和精度,模擬了電磁波垂直入射到任意磁偏角磁化等離子體平板的反射/透射系數(shù)幅度和相位. 入射電磁波采用微分高斯脈沖,沿+z軸傳播3.75 ns,傳播模型如圖1. 計(jì)算空間占800個(gè)網(wǎng)格,等離子體占第400~600個(gè)網(wǎng)格,均勻分布,厚1.5 cm,網(wǎng)格兩端采用一階Mur吸收邊界條件[16],其余為自由空間. 根據(jù)Courant穩(wěn)定性條件,空間迭代步長Δz=75 μm,迭代時(shí)間Δt=0.125 ps. 等離子體參數(shù)設(shè)置:碰撞頻率vc=20 GHz,電子回旋頻率ωce=16 GHz,等離子體頻率ωp=28.7 GHz.

        圖1 電磁波沿任意磁偏角等離子體層的傳播模型

        記錄第399和第601個(gè)網(wǎng)格的電場值,進(jìn)行離散傅里葉變換,可從式(8)得到特征波頻域內(nèi)的反射、透射系數(shù):

        R(ω)=[Erx(ω)+κ·Ery(ω)]/Eix(ω),

        (8)

        T(ω)=[Etx(ω)+κ·Ety(ω)]/Eix(ω).

        (9)

        式中:

        Erx(ω)=DFT[Etx(t)-Eix(t)];

        Ery(ω)=DFT[Ety(t)-Eix(t)];

        Etx(ω)=DFT[Etx(t)];

        Ety(ω)=DFT[Ety(t)];

        (10)

        電磁波傳播方向與外磁場垂直,“橫向傳播”(θ=π/2),特征波橫向復(fù)數(shù)比取值分兩種情況:κ→0,E矢量或它在xy平面內(nèi)的分量所描繪的橢圓退化成直線(y方向線極化)的O波,電磁波傳播不受磁場影響;κ→∞,Ey=0,E矢量為位于xz平面內(nèi)(橢圓極化)的X波. 圖2是尋常波O波和非尋常波X波的電場矢量E極化示意圖,O波電場矢量和磁場平行,X波電場矢量和磁場垂直.

        (a)尋常波O波 (b)非尋常波X波圖2 尋常波和非尋常波電場矢量極化圖

        圖3是磁偏角θ=π/2,O波的反射/透射系數(shù)幅度和相位的結(jié)果. O波反射/透射系數(shù)幅度和相位的數(shù)值計(jì)算結(jié)果與解析解吻合. 當(dāng)入射電磁波頻率ω>28.7 GHz時(shí),透射系數(shù)幅度大于反射系數(shù)幅度,電磁波易傳播,這與O波只有一個(gè)傳播帶,滿足電磁波頻率大于等離子體頻率ω>ωp相吻合.

        圖4是θ=π/2,X波的反射/透射系數(shù)幅度和相位的結(jié)果. X波反射/透射系數(shù)幅度和相位的數(shù)值計(jì)算結(jié)果也與解析解吻合. 當(dāng)入射電磁波頻率分別滿足0<ω<21.8 GHz和32.8 GHz<ω<37.7 GHz時(shí),這兩個(gè)截止區(qū)反射系數(shù)幅度較大,電磁波傳播困難;入射電磁波頻率分別滿足21.8 GHz<ω<32.8 GHz和ω>37.7 GHz時(shí),透射系數(shù)幅度較大,電磁波易傳播. 原因在于X波有兩個(gè)截止帶0<ω<ωL和ωuh<ω<ωR,兩個(gè)傳播帶ωL<ω<ωuh和ω>ωR[18],ωR為右旋圓極化X波截止頻率,ωL為左旋圓極化X波截止頻率,ωuh為高混雜頻率,分別滿足式(11)~(13). 數(shù)值計(jì)算得到的波傳播和截止帶與理論計(jì)算一致.

        (11)

        (12)

        (13)

        圖3 O波反射/透射系數(shù)幅度和相位

        圖4 X波反射/透射系數(shù)幅度和相位

        2.1.2 法拉第旋轉(zhuǎn)效應(yīng)

        線極化電磁波穿過磁化等離子體,可以分解成兩個(gè)不同相速度的圓極化波,它們的色散關(guān)系可以表示成[19]:

        (14)

        (15)

        式中:kL,kR分別為左、右旋波傳播常數(shù);c為光速.假設(shè)入射電磁波為x方向極化的單頻線極化波,則透過厚度為d的等離子體后的線極化波與x方向的傾斜角Ω即法拉第旋轉(zhuǎn)角可表示為

        (16)

        為了說明法拉第旋轉(zhuǎn)效應(yīng),用任意磁偏角JEC-FDTD仿真求解θ=0°,入射電磁波頻率10 GHz的x方向線極化波,等離子體參數(shù)ωp=5 GHz,ωce=16 GHz,vc=5 GHz,不同等離子體板厚度(0,100,200,300,400個(gè)網(wǎng)格)透射場,得到電場軌跡如圖5所示. 結(jié)果表明:隨著等離子體厚度的增加,法拉第旋轉(zhuǎn)角增加,電場強(qiáng)度軌跡的橢圓度也越大.這是由于左旋和右旋波的相速度不同以及電磁波的碰撞衰減造成的,這一現(xiàn)象與理論分析吻合.

        圖5 不同厚度磁化等離子體透射場幅度軌跡

        2.2 任意磁偏角傳播特性

        前面對于電磁波在磁化等離子體中傳播(θ=π/2)和法拉第旋轉(zhuǎn)特性進(jìn)行了分析(θ=0°),本節(jié)給出任意磁偏角情況下,電磁波的傳播特性.吸收率的定義滿足式(17),R、T定義與前面相同[20].

        A=1-R2-T2.

        (17)

        圖6給出了電磁波沿任意磁偏角(θ=0°,30°,60°, 90°)等離子體傳播時(shí)兩種波的反射/透射系數(shù)幅度和吸收率的大小. 圖6(a)是I波(式(10)取“-”)在任意磁偏角情況下的反射/透射系數(shù)幅度和吸收率的大小,當(dāng)θ=0°(κ=-j),為左旋圓極化(Left Hand Circularly Polarized,LCP)波. 圖6(b)是II波(式(10)取“+”)在任意磁偏角情況下的反射/透射系數(shù)幅度和吸收率大小,θ=0°(κ=+j),為右旋圓極化(Right Hand Circularly Polarized,RCP)波. LCP波有一個(gè)截止帶, 一個(gè)傳播帶,當(dāng)ω>ωL才能傳播,而RCP波則有兩個(gè)傳播帶0<ω<ωce,ω>ωR和一個(gè)截止帶ωce<ω<ωR.ωL,ωR分別是LCP波和RCP波的截止頻率[18],與前面X波左右旋圓極化波截止頻率定義相同,其他各參數(shù)設(shè)置也與前面算例相同. 結(jié)果分析類似于前面X波的傳播,不同的是RCP波只有一個(gè)截止帶而X波卻有兩個(gè).

        (a) I波(式(10)取“-”)

        (b) II波(式(10)取“+”)圖6 任意磁偏角I、II波的反射/透射系數(shù)幅度和吸收率大小曲線

        圖6結(jié)果表明:入射電磁波頻率大約在小于24 GHz時(shí),隨著磁偏角的增大I波衰減增加,II波則相反,等離子體吸收減弱;I波和II波在等離子頻率 (28.7 GHz)附近,反射/透射系數(shù)幅度有明顯的突變, 說明等離子頻率是一個(gè)反射點(diǎn); 高頻區(qū)波的衰減都較小,更適合兩種波的傳播;隨著θ角增加波的吸收帶寬都變大.

        2.3 等離子體共振吸收

        2.3.1 電子回旋共振吸收

        高頻入射電磁波傳播方向與背景磁場平行ki‖B(θ=0°),色散關(guān)系可分別表示[18]為:

        (18)

        (19)

        圖7給出了RCP/LCP波吸收率與等離子體頻率和碰撞頻率的關(guān)系.參數(shù)設(shè)置分別為(a)ωce=32 GHz,ωp=10, 28.7, 50 GHz,vc=20 GHz;(b)ωce=32 GHz,ωp=28.7 GHz,vc=10, 20, 30 GHz. RCP波吸收率只在電子回旋頻率附近出現(xiàn)峰值,產(chǎn)生共振吸收;隨著等離子體頻率的增加,共振吸收峰值左移,吸收帶寬增加. 隨著碰撞頻率的增加,共振峰值逐漸變寬. LCP波沒有發(fā)生共振吸收,這是因?yàn)樗惴ê雎粤穗x子運(yùn)動(dòng).

        (a) 電子回旋共振吸收與等離子體頻率關(guān)系

        (b) 電子回旋共振吸收與碰撞頻率關(guān)系圖7 RCP/LCP波電子回旋共振吸收與等離子體頻率、碰撞頻率關(guān)系

        2.3.2 高頻混雜共振吸收

        電磁波傳播方向與背景磁場垂直(θ=π/2),磁化等離子體中的X波無碰撞情況下滿足色散關(guān)系[18]

        (20)

        式中:ωL,ωR,ωuh定義與前面相同;“+”為右旋圓極化波.由式(11)~(13)顯然能夠得到ωL<ωuh<ωR,X波中右旋圓極化波kE為無窮純虛數(shù),使該波傳播很短距離振幅就迅速衰減到很小,產(chǎn)生高頻混雜共振吸收.

        圖8給出了X波/O波吸收率隨電子回旋頻率,等離子體碰撞頻率的變化關(guān)系.參數(shù)設(shè)置為(a)ωce=16, 32, 48 GHz,ωp=28.7 GHz,vc=20 GHz,可得高頻混雜波頻率分別為ωuh=32.9, 43, 55.9 GHz;(b)ωp=28.7 GHz,ωce=32 GHz,vc=10, 20, 30 GHz. 從圖8可以得到:X波吸收率在高頻混雜頻率處出現(xiàn)峰值;而且碰撞頻率越大,吸收峰值帶寬越大. 由于等離子體朗繆爾共振,使O波吸收率在等離子體頻率附近突變,產(chǎn)生吸收峰.

        (a) 高頻混雜共振吸收與電子回旋頻率關(guān)系

        (b) 高頻混雜共振吸收與碰撞頻率關(guān)系圖8 X/O波高頻混雜共振吸收與電子回旋頻率、碰撞頻率關(guān)系

        2.4 電離層加熱等離子體高頻混雜共振吸收

        為了模擬實(shí)際電離層加熱實(shí)驗(yàn),計(jì)算了空間400個(gè)網(wǎng)格,等離子體占150~350網(wǎng)格,60 m厚,迭代步長0.3 m,迭代時(shí)間500 ps.依據(jù)高緯度地區(qū)高頻極光計(jì)劃HAARP電離層加熱典型參數(shù)[2,22-23],入射電磁波選擇頻率1~10 MHz的微分高斯脈沖,電子回旋頻率ωce=1.4 MHz, 等離子體頻率ωp=2.9, 4.3, 5.7, 7.1 MHz,考慮碰撞頻率vc=50 kHz.

        圖9給出了入射波頻率在nωce+0.1(n=2~5)處的X/O波高頻混雜共振吸收結(jié)果, 圖中黑色虛線從左到右所在頻率分別為nωce,ωp,ω0,ωuh. 仿真結(jié)果表明ωp<ω0<ωuh時(shí),共振吸收最強(qiáng),與實(shí)驗(yàn)觀測到激發(fā)人工電離層和受激上行拓展輻射譜的頻率關(guān)系保持一致[23-24].

        圖9 不同電子回旋諧波加熱(n=2~5)時(shí)X/O波的高頻混雜共振吸收現(xiàn)象

        3 結(jié) 論

        本文給出任意磁偏角JEC-FDTD算法的推導(dǎo),計(jì)算一維磁化等離子體反射/透射系數(shù)幅度和相位以及法拉第旋轉(zhuǎn)效應(yīng),通過與解析解比較驗(yàn)證了算法的有效性、精確性. 研究了等離子體共振吸收特性,結(jié)果表明:隨著磁偏角增加,電磁波在等離子頻率處的突變特性增強(qiáng);RCP波在電子回旋諧振頻率處發(fā)生共振,吸收系數(shù)最大;非尋常波X波在高頻混雜頻率處發(fā)生共振,吸收最大. 背景磁場和電子密度分布會(huì)影響共振吸收峰位置,增加等離子碰撞頻率可以增大電磁波的吸收帶寬. 因此,選擇合適的背景磁場、電子密度和碰撞頻率可以使電磁波達(dá)到更好的吸收.重點(diǎn)研究了電離層加熱高頻混雜共振吸收的定量計(jì)算,結(jié)果與實(shí)驗(yàn)觀測到的現(xiàn)象一致,在加熱頻率滿足高頻混雜共振吸收條件時(shí),電磁波與等離子體的非線性作用使其能量更多地被轉(zhuǎn)換成等離子體能量,提高了電離層加熱效率.

        實(shí)際電離層環(huán)境不僅是任意磁偏角, 等離子體的電子密度等參數(shù)還是隨時(shí)間不斷變化的,時(shí)變非均勻的等離子體變化也會(huì)影響電磁波傳播和吸收,電離層加熱的效率和實(shí)際試驗(yàn)場景也是密切相關(guān)的[25-27]. 時(shí)變等離子體需要采用更為復(fù)雜的FDTD迭代, 把電子密度隨時(shí)間變化關(guān)系考慮到迭代中去,這是今后進(jìn)一步要研究的內(nèi)容.

        [1] FU H Y, SCALES W A, BERNHARDT P A, et al. Stimulated Brillouin scattering during electron gyro-harmonic heating at EISCAT[J]. Annales geophysicae, 2015, 33:983-990.

        [2] FU H Y, SCALES W A, BERNHARDT P A, et al. Stimulated Brillouin scatter and stimulated ion Bernstein scatter during electron gyroharmonic heating experiments[J]. Radio science, 2016, 48(5):607-616.

        [3] LAROUSSI M, ROTH J R. Numerical calculation of the reflection, absorption, and transmission of microwaves by a nonuniform plasma slab[J]. IEEE transactions on plasma science, 2002, 21(4):366-372.

        [4] HU B J, WEI G, LAI S L. SMM analysis of reflection, absorption, and transmission from no uniform magnetized plasma slab[J]. IEEE transactions on plasma science, 1999, 27(4):1131-1136.

        [5] 柳文, 焦培南, 王俊江, 等. 利用射線追蹤研究電離層擾動(dòng)[J].地球物理學(xué)報(bào), 2005, 48(3): 465-470.

        LIU W, JIAO P N, WANG J J, et al. A study of the disturbance in the ionosphere using ray tracing technology[J]. Chinese journal of geophysics, 2005, 48(3): 465-470.(in Chinese)

        [6] CROFT T A, HOOGANSIAN H. Exact ray calculations in a quasi-parabolic ionosphere with no magnetic field[J]. Radio science, 1968, 3(1): 69-74.

        [7] SULLIVAN D M. Frequency-dependent FDTD methods using Z transform [J]. IEEE transactions on antennas & propagation, 1992, 40(10):1223-1230.

        [8] SHIBAYYAMA J, ANDO R, NOMURA A, et al. Simple trapezoidal recursive convolution technique for the frequency-dependent FDTD analysis of a Drude-Lorentz model[J]. IEEE photonics technology letters, 2009, 21(2):100-102.

        [9] LIU S, CHEN S, HE Z, et al. Study on the polarization properties of electromagnetic waves with arbitrary magnetic declination in magnetized plasma[J]. Waves in random and complex media, 2015, 25(3):1-12.

        [10] ATTIYA A M, ABDULLAH H H. Shift-operator finite difference time domain: an efficient unified approach for simulating wave propagation in different dispersive media[J]. IEEE antennas and propagation, 2010:1-4.

        [11] WANG F, WEI B, GE D B. A method for FDTD modeling of wave propagation in magnetized plasma[C]//International Conference on Consumer Electronics, Communications and Networks. IEEE, 2011:4659-4662.

        [12] MA L X, ZHANG H, ZHENG H X, et al. Shift-operator FDTD method for anisotropic plasma in KDB coordinates system[J]. Progress in electromagnetics research M, 2010, 12(11): 51-65.

        [13]CHEN Q, KATSURAI M, AOYAGI P H. An FDTD formulation for dispersive media using a current density[J]. IEEE transactions on antennas and propagation, 1998, 46(11):1739-1746.

        [14] XU L J, YUAN N C. FDTD formulations for scattering from 3-D anisotropic magnetized plasma objects[J]. IEEE antennas and wireless propagation letters, 2006, 5(1):335-338.

        [15] QIAN Z H, CHEN R S. FDTD Analysis of magnetized plasma with arbitrary magnetic declination[J]. Journal of infrared, millimeter, and terahertz waves, 2007, 28(2): 157-167.

        [16] TAFLOVE A. Advances in computational electrodynamics: the finite-difference time-domain method[M]. Boston: Artech House, 1998.

        [17] GINZBURG V L, TAYLER R J. Physics. (book reviews: the propagation of electromagnetic waves in plasmas)[J]. Science, 1965: 147.

        [18] CHEN F F. Introduction to plasma physics and controlled fusion[J]. Physics today, 1985, 38(5): 87-88.

        [19] BOOKER H G. Cold plasma waves[M]. Hingham: Kluwer, 1984.

        [20] WANG M Y. Propagation properties of terahertz waves in a time-varying dusty plasma slab using FDTD[J]. IEEE transactions on plasma science, 2015, 43(12): 4182-4186.

        [21] SHARMA A S, ELIASSON B, MILIKH G M, et al. Low-frequency waves in HF heating of the ionosphere[M]//Low-frequency waves in space plasmas. Hoboken: John Wiley & Sons, Inc., 2016.

        [22] GUREVICH A V. Nonlinear effects in the ionosphere[J]. Physics-uspekhi, 2007, 50(11): 1145-1177.

        [23] SERGEEV E, GRACH S, SHINDIN A, et al. Artificial ionospheric layers during pump frequency stepping near the 4th gyroharmonic at HAARP[J]. Physical review letters, 2013, 110(6): 065002.

        [24] LEYSER T B, THIDé B, DERBLOM H, et al. Dependence of stimulated electromagnetic emission on the ionosphere and pump wave[J]. Journal of geophysical research space physics, 1990, 95(A10):17233-17244.

        [25] CHENG M S, XU B, WU Z S, et al. Observation of VHF incoherent scatter spectra disturbed by HF heating[J]. Journal of atmospheric and solar-terrestrial physics, 2011, 105/106: 245-252.

        [26] XU B, CHENG M S, XU Z W, et al. The analysis of an ionospheric heating experiment in polar region[C]//2014 XXXIth URSI General Assembly and Scientific Symposium. Beijing: IEEE, 2014:1-4.

        [27] 徐彬,王占閣,許正文,等.極區(qū)冬季電離層加熱實(shí)驗(yàn)研究(三)——低電離層分析[J].地球物理學(xué)報(bào),2010,53(6):1263-1268.

        XU B, WANG Z G, XU Z W, et al. Observations of the heating experiments in the polar winter ionosphere III-analysis in the low region[J]. Chinese journal of geophysics, 2010, 53(6): 1263-1268.(in Chinese)

        日本激情视频一区在线观看| 久久综合精品国产一区二区三区无码| 成年女人毛片免费观看97| 精品福利一区| 精品人妻av区二区三区| 午夜视频在线观看一区二区小| 国产乱国产乱老熟300部视频| 亚欧AV无码乱码在线观看性色 | 精品人妻伦九区久久aaa片| 自拍偷自拍亚洲精品播放| 成人国产一区二区三区精品不卡 | 亚洲第一网站免费视频| 久久精品国产亚洲av电影网| 久久中文字幕乱码免费| 一本久道久久综合狠狠操| 激情在线一区二区三区视频| 成l人在线观看线路1| 在线一区不卡网址观看| 免费无遮挡无码永久视频| 亚洲熟女少妇精品综合| 少妇人妻在线无码天堂视频网| 久久精品日韩av无码| 男女上床视频在线观看| 亚洲国产一区二区三区精品| 小sao货水好多真紧h无码视频| 专区国产精品第一页| 美女草逼视频免费播放| 久久无码高潮喷水抽搐| 中文字幕无码av激情不卡| 亚洲av在线播放观看| 久久人妻少妇嫩草av蜜桃| 国产专区一线二线三线码 | 亚洲一区二区三区在线激情| 人妻少妇-嫩草影院| 色狠狠一区二区三区香蕉| 欧美巨大xxxx做受中文字幕| 国产精品白浆一区二区免费看| 国产精品第一国产精品| 手机看片国产日韩| 久亚洲一线产区二线产区三线麻豆 | 九九日本黄色精品视频|