999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

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

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

張景 付海洋

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

引 言

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

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

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

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

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

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

(1)

(2)

(3)

(4)

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

(5)

(6)

(7)

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

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

2.1 算例驗證

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

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

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

記錄第399和第601個網(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ù)值計算結(jié)果與解析解吻合. 當(dāng)入射電磁波頻率ω>28.7 GHz時,透射系數(shù)幅度大于反射系數(shù)幅度,電磁波易傳播,這與O波只有一個傳播帶,滿足電磁波頻率大于等離子體頻率ω>ωp相吻合.

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

(11)

(12)

(13)

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

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

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

線極化電磁波穿過磁化等離子體,可以分解成兩個不同相速度的圓極化波,它們的色散關(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個網(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ù)幅度和吸收率的大小. 圖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波有一個截止帶, 一個傳播帶,當(dāng)ω>ωL才能傳播,而RCP波則有兩個傳播帶0<ω<ωce,ω>ωR和一個截止帶ωce<ω<ωR.ωL,ωR分別是LCP波和RCP波的截止頻率[18],與前面X波左右旋圓極化波截止頻率定義相同,其他各參數(shù)設(shè)置也與前面算例相同. 結(jié)果分析類似于前面X波的傳播,不同的是RCP波只有一個截止帶而X波卻有兩個.

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

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

圖6結(jié)果表明:入射電磁波頻率大約在小于24 GHz時,隨著磁偏角的增大I波衰減增加,II波則相反,等離子體吸收減弱;I波和II波在等離子頻率 (28.7 GHz)附近,反射/透射系數(shù)幅度有明顯的突變, 說明等離子頻率是一個反射點; 高頻區(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ā)生共振吸收,這是因為算法忽略了離子運動.

(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 電離層加熱等離子體高頻混雜共振吸收

為了模擬實際電離層加熱實驗,計算了空間400個網(wǎng)格,等離子體占150~350網(wǎng)格,60 m厚,迭代步長0.3 m,迭代時間500 ps.依據(jù)高緯度地區(qū)高頻極光計劃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時,共振吸收最強(qiáng),與實驗觀測到激發(fā)人工電離層和受激上行拓展輻射譜的頻率關(guān)系保持一致[23-24].

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

3 結(jié) 論

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

實際電離層環(huán)境不僅是任意磁偏角, 等離子體的電子密度等參數(shù)還是隨時間不斷變化的,時變非均勻的等離子體變化也會影響電磁波傳播和吸收,電離層加熱的效率和實際試驗場景也是密切相關(guān)的[25-27]. 時變等離子體需要采用更為復(fù)雜的FDTD迭代, 把電子密度隨時間變化關(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] 柳文, 焦培南, 王俊江, 等. 利用射線追蹤研究電離層擾動[J].地球物理學(xué)報, 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ū)冬季電離層加熱實驗研究(三)——低電離層分析[J].地球物理學(xué)報,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)

主站蜘蛛池模板: 九九九久久国产精品| 国产第一页免费浮力影院| 日韩在线观看网站| 亚洲欧美日韩精品专区| 国产欧美中文字幕| 国产美女91视频| 亚洲精品天堂自在久久77| 国产高清不卡视频| JIZZ亚洲国产| 国产日韩精品一区在线不卡| 人妻精品久久久无码区色视| 久久综合九色综合97婷婷| 伊人久久影视| 亚洲中文字幕在线观看| 亚洲精品国产首次亮相| h视频在线播放| 毛片手机在线看| 午夜爽爽视频| AV老司机AV天堂| 人妻少妇乱子伦精品无码专区毛片| 国产99久久亚洲综合精品西瓜tv| 2018日日摸夜夜添狠狠躁| 亚洲精品麻豆| 欧美亚洲一区二区三区导航| 波多野吉衣一区二区三区av| 国内a级毛片| 日本午夜网站| 久久亚洲美女精品国产精品| 久久亚洲国产一区二区| 久久国产精品无码hdav| 日本高清有码人妻| 久久精品无码国产一区二区三区| 国产精品欧美日本韩免费一区二区三区不卡| 久久永久免费人妻精品| WWW丫丫国产成人精品| 亚洲最大福利网站| 全部免费毛片免费播放| 国产成人av一区二区三区| 国产福利免费观看| 欧美激情第一区| 国产黄视频网站| 国产精品香蕉| 毛片一级在线| av午夜福利一片免费看| 婷婷五月在线视频| 国产裸舞福利在线视频合集| 亚洲无码91视频| 99精品国产高清一区二区| 国产天天射| 1024你懂的国产精品| 天天色天天综合网| 秋霞午夜国产精品成人片| 91精品久久久无码中文字幕vr| 试看120秒男女啪啪免费| 欧美成人精品欧美一级乱黄| 亚洲国产日韩在线成人蜜芽| 国产免费福利网站| 亚洲成在线观看| 超碰91免费人妻| 国产精品亚洲αv天堂无码| 中文字幕永久视频| 最新国语自产精品视频在| 97视频免费看| 欧美午夜在线播放| 日韩欧美网址| 亚洲AV无码不卡无码| 国产91高跟丝袜| 夜夜拍夜夜爽| 欧美日韩精品一区二区视频| 亚洲成综合人影院在院播放| 丁香五月婷婷激情基地| 蝴蝶伊人久久中文娱乐网| 成·人免费午夜无码视频在线观看| 久久99国产综合精品1| 久久窝窝国产精品午夜看片| 国产va免费精品| 精品国产免费观看一区| 99热这里都是国产精品| 老司国产精品视频91| 亚洲VA中文字幕| 亚洲天堂首页| 99伊人精品|