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

氙氣中輻射激波的發(fā)光特性

2021-05-07 06:08:14趙多李守先安建祝2吳勇吳澤清李瓊王芳孟廣為
物理學報 2021年7期
關鍵詞:實驗模型

趙多 李守先 安建祝2)? 吳勇 吳澤清李瓊 王芳 孟廣為

1) (北京應用物理與計算數學研究所, 北京 100094)

2) (北京大學應用物理與技術研究中心, 北京 100871)

本文模擬研究了氙氣中X 射線加熱產生輻射激波的發(fā)光特性.輻射激波采用Zinn 模型計算, 并在模型中輸入氙氣的輻射不透明度和狀態(tài)方程參數.研究發(fā)現, 輻射激波伴隨著豐富的光學演化過程, 激波對外輻射強度表現出兩個明顯的亮度峰值和一個亮度極小值, 輻射光譜也經常偏離黑體輻射光譜.對氙氣不同位置光學特性的分析可知, 激波和內部高溫區(qū)的輻射吸收系數的動態(tài)演化導致了輻射激波的發(fā)光位置和輻射強度的變化.

1 引 言

輻射激波一般是非常強的激波, 波陣面處介質強壓縮導致的輻射效應已經可以影響激波的宏觀性質, 需要采用輻射流體力學理論描述這一現象.輻射激波在天體物理和慣性約束聚變中都受到廣泛關注.在超新星爆發(fā)[1]、天體噴流[2]、吸積盤[3]等活動中, 輻射對激波的能量損失和流體不穩(wěn)定性的發(fā)展有重要影響.在間接驅動慣性約束聚變中, 靶丸受X 射線照射產生的強激波也伴隨輻射效應,導致能量的損失和激波結構的改變[4].

近年來, 利用大型激光裝置創(chuàng)造的高能量密度條件可以產生高強度的輻射激波[5?7], 開展相關天體物理現象的實驗模擬研究[8].實驗中一般通過激光輻照固體靶產生高溫等離子體, 靶物質壓縮或加熱實驗氣體產生輻射激波[9,10].不同靶物質與實驗氣體的作用方式不同, 塑料靶主要通過靶物質激波推動和壓縮氣體產生輻射激波; 而高Z 金屬靶輻射的X 射線可以快速加熱周圍的氣體, 演化出輻射激波.實際應用中經常設計混合材料的固體靶,以產生特定的輻射激波[11].目前開展的實驗中已經對輻射激波的形態(tài)結構[12,13]、輻射光譜[14,15]、激波參數[16?18]和不穩(wěn)定性發(fā)展[6,19]等特征進行了廣泛的研究.利用輻射激波產生的高溫高壓條件也可以研究極端狀態(tài)下的物質性質[20], 并獲取相應的狀態(tài)方程[21]和輻射不透明度[22]數據.

輻射激波的理論研究也一直在發(fā)展中.Sew和Guess[23]發(fā)現輻射效應將導致激波陣面的展寬.Zel’dovich 和Raizer[24]指出激波發(fā)出的輻射將加熱上游的介質形成波前結構, 并導致激波陣面后形成一個溫度尖峰.Drake[25]使用半解析方法發(fā)現光學厚介質中輻射激波的波前包括擴散區(qū)和透射區(qū),且擴散區(qū)的特征不完全依賴于激波下游的性質.McClarren 等[26]指出當激波上游為光學薄介質時,波前不會形成, 激波發(fā)出的輻射將離開系統(tǒng), 導致激波后的密度要高于光學厚的情況.Mabey 等[27]發(fā)展了光學厚介質中輻射激波下游溫度的計算方法, 為實驗中激波溫度的評估提供依據.

本文使用模擬方法研究輻射激波發(fā)光特性的動態(tài)演化及其形成機制.目前的實驗研究難以測量輻射激波內部的光學特征, 而理論研究多局限于特定條件下的激波性質, 很少關注激波的演化過程,已發(fā)表文獻中對輻射激波光學特征的演化多是一些定性描述, 缺乏定量、清晰的研究結果.本研究關注X 射線加熱氙氣形成的輻射激波, 在這一條件下激波的演化過程和發(fā)光現象非常豐富, 研究結果有助于理解輻射激波的物理過程,并為相關實驗設計提供理論支撐.本文第2 節(jié)分別介紹計算模型(2.1), 計算程序中輸入的氙氣參數(2.2), 以及程序和參數的校驗方法(2.3); 第3 節(jié)是本文的研究結果和討論; 第4 節(jié)是全文總結.

2 數值模擬方法與程序校驗

2.1 計算模型

本文對輻射激波的模擬使用Zinn[28]提出的一維輻射流體力學模型(本文稱為Zinn 模型).該模型的氣體動力學和輻射輸運過程均在球對稱的拉格朗日坐標下求解.模型假設介質處于局部熱平衡狀態(tài), 并采用多群近似求解輻射輸運過程.這一方法在輻射能流的計算中充分利用了計算網格的球對稱特征, 并對高溫區(qū)和周圍介質的輻射能流區(qū)分處理, 大大加快了計算速度.

Zinn 模型采用有限差分格式, 將空間劃分為一系列同心球層, 并假設每一球層的密度、溫度等力學量是均勻分布的.輻射輸運過程求解中, 由P1點傳輸到P2點輻射的譜強度可以表示為(由忽略散射項的積分形式輻射輸運方程求得)

其中I?為輻射的譜強度,?表示輻射頻率, ?s為P1和P2間的距離,μ′a為考慮受激發(fā)射的輻射吸收系數,B?為平衡輻射的譜強度.

輻射的譜能流Fυ定義在球層邊界, 同一位置向外和向內的譜能流可分別表示為:

其中θ為半徑方向與特定輻射方向的夾角.將I?的表達式代入, 并根據球對稱幾何關系可以得到任一球層通過內邊界向內的輻射譜能流F1?和通過外邊界向外的輻射譜能流F2+(F1?和F2+均略去下標?)的表達式.由于計算過程比較繁瑣, 此處略去具體的表達式和推導過程, 這些都可在文獻[28]中找到.

本模型的一個重要假設是空間中每一點的輻射譜強度都可以用雙值階梯函數來表示, 即來自核心高溫區(qū)的輻射取大的強度值Ia, 核心區(qū)以外的輻射取較小強度值Ib.將這兩個輻射譜強度值代入球層內外邊界輻射譜能流的表達式, 聯(lián)立消去輻射譜強度變量便可得到F1?與F2+的遞推關系.其中F1?僅與其外層向內的輻射譜能流F2?有關, 可以從最外層(輻射譜能流一般為零)逐層向內求解.將計算得到的一系列向內的輻射譜能流代入, 并從核心(半徑和輻射譜能流均為零)逐層向外求解F2+,這樣就得到了系統(tǒng)對外的輻射譜能流.當兩個相鄰的球層均為光學厚介質(光學厚度大于2)時, 上述輻射輸運過程可以簡化為擴散近似來求解.

本模型計算每一時間步輻射輸運的同時, 進行一步或多步流體計算.流體求解過程中, 首先根據密度、內能和狀態(tài)方程求解每一球層的壓強Pi和粘性壓Qi, 然后根據這些壓強值求解每一層在流體時間步 ?tH中產生的加速度, 通過加速度可以計算該層在 ?tH后的位移、速度和內能, 進一步得出其他流體參量.

2.2 氣體參數

在模型計算過程中需要輸入氣體的狀態(tài)方程與輻射不透明度參數, 這些參數的準確性對于得到可靠的計算結果非常重要.本文使用的氙氣輻射不透明度數據采用DCA/UTA 不透明度程序計算[29].該程序中原子的能級結構采用細致組態(tài)近似(detailed configuration accounting, DCA), 在計算譜線躍遷時, 使用UTA(unresolved transition array)統(tǒng)計模型近似考慮細致能級的譜線結構[30,31].這套程序計算的不透明度結果與其他理論和實驗結果進行了大量比較[29,32], 證實其計算結果是可信的,且具有較高的精度.本文計算得到兩個密度下的氙氣不透明度數據如圖1 所示.

氙氣狀態(tài)方程數據的計算結合化學模型[33]與平均原子模型, 其中平均原子模型的原子結構用哈特里-福克-斯萊特(HFS)自洽場方法計算[34].計算中, 當溫度低于2.6 eV 時采用化學模型, 溫度高于8.6 eV 時采用平均原子模型, 兩個溫度之間的區(qū)域采用插值數據.化學模型基于經驗勢和液體變分微擾理論計算了氙分子間相互作用能, 考慮了四級電離, 并基于正則分布和NIST 公開數據庫[35]提供的原子激發(fā)能級計算了離子組分內部的電子熱激發(fā)能.計算得到的狀態(tài)方程參數如圖2 所示.

2.3 模型校驗

為了驗證所用程序和氙氣物性參數的可靠性,本文模擬了Vinci 等[17]和Koenig 等[18]對氙氣中輻射激波的實驗測量結果.這些實驗中激波產生的原理是: 激光照射到靶材上, 產生燒蝕激波, 靶材的設計避免了X 射線加熱, 這樣燒蝕激波作為活塞, 推動氙氣產生輻射激波.實驗中分別測量了0.1 和0.2 個標準大氣壓(atm)下氙氣中輻射激波的平均速度和溫度.

Zinn 模型對0.1 atm 氙氣的計算中, 輻射激波的溫度在1 ns 時達到最大值, 約為20 eV (由380—460 nm 波段的輻射反推, 其他可見光波段的計算結果差異不大), 在實驗測量中這一數值約為18 eV (文獻[17]的圖3), 模擬計算2—4 ns 的激波平均速度為63 km/s, 與實驗測量值65 km/s也很接近.在0.2 atm 的計算中, 輻射激波的最高溫度和平均速度分別為13 eV 和38 km/s, 與實驗中的兩次測量值9 eV, 45 km/s 和11 eV, 41 km/s(文獻[18]的表2)相差不多.以上對比表明, 本文使用的模型和氣體參數在實驗涉及的氙氣物性參數范圍內是較為可靠的.

圖1 不同氣體密度時氙氣的輻射不透明度隨光子能量的變化 (a) ρ = 5.33 × 10–7 g/cm3; (b) ρ = 5.33 × 10–3 g/cm3.不同曲線類型表示不同的溫度 (1 eV = 11610 K)Fig.1.Opacity data of xenon at the density of 5.33 × 10–7 g/cm3 (a) and 5.33 × 10–3 g/cm3 (b), lines with different style represent different temperature.

圖2 氙氣的狀態(tài)方程數據 (a)溫度(T, 單位為eV)與比內能(E, 單位為erg/g (1 erg/g = 10–7 J/g))之比隨E 的變化; (b)壓強(P, 單位dyn/cm2 (1 dyn/ cm2 = 10–5 N/ cm2))與E 和密度( ρ , 單位g/cm3)乘積之比隨E 的變化.參數使用比值形式是為方便程序計算Fig.2.Equation-of-state data of xenon: (a) Variation of the ratio between temperature (T, in eV) and specific inertial energy (E, in erg/g) with E; (b) variation of ratio between pressure (P, in dyn/cm2) and the multiplication of E and density ( ρ , in g/cm3) with E.The parameters are shown by ratios for computation convenience.

圖3 不同時刻氙氣中的密度(a)和溫度(b)的位置分布Fig.3.Variation of density (a) and temperature (b) with position at different times in xenon.

Zinn 模型不包含激光與靶物質的相互作用過程, 因此在計算中對實驗情形做了簡化.激波通過最內層球面的高速膨脹產生(相當于燒蝕激波的推動), 形成穩(wěn)定的激波后, 計算激波各參量并與實驗情況作對比.在實驗中, 氙氣激波既有在燒蝕激波推動下向前的運動, 又有橫向的擴張.本文模型僅描述一維球對稱擴張的情形, 并選取最內層驅動球面的半徑為0.15 cm, 與實驗情況類比相當于某個立體角所截取的波陣面的運動和擴張過程.模擬中激波的橫向尺度在7 ns 內擴張1.33 倍, 實驗中的橫向擴張是1.67 倍.原文中也曾使用一維程序進行模擬, 即使完全不考慮橫向擴張過程, 仍然可以模擬激波的發(fā)光特征[17,18], 因此就一維模擬而言,本文計算更接近實驗中輻射激波的擴張過程.模型計算的多個參量都與實驗測量值接近也證明了模型對激波物理過程描述的合理性.

3 結果與討論

本文研究的輻射激波由X 射線加熱氙氣產生.模擬設置中, 激光能量以內能的形式加載到核心球內, 核心球的半徑取0.03 cm, 能量加載使用高斯函數, 總能量為250 J, 持續(xù)時間為1 ns.在這一條件下, 核心球的溫度可以被加熱到100 eV 以上.本研究不涉及激光與靶物質的耦合過程, 但通過能量加載函數來代替靶物質X 射線向氣體中傳輸能量的過程.對于不同初始條件的計算表明, 只要總能量不變且核心球能夠被加熱到100 eV 以上, 激光能量加載函數和核心球半徑的小幅調整, 對輻射激波演化過程的影響不大, 表明本文的計算結果是穩(wěn)定的.核心球和周圍氙氣密度均取2.58 × 10–2g/cm3, 對應常溫時的壓力約為4.5 atm.以上條件參考神光II 裝置的激光參數[5].

為了研究能量輸入后氙氣中激波的演化過程,圖3 給出了幾個典型時刻激波的密度(圖3(a))和溫度(圖3(b))分布.早期(< 3 ns)的氣體密度擾動很小, 激波還未形成, 但中心溫度高達幾十電子伏(eV), 形成一個高溫區(qū).此時高溫區(qū)的擴張并未伴隨明顯的流體效應, 主要通過輻射加熱周圍介質, 形成高速擴張的熱波[36].一段時間后(20 ns),流體擾動趕上熱波陣面并以激波的形式進入周圍的氙氣中, 此時激波陣面與高溫區(qū)表面重合, 高溫區(qū)的密度明顯降低.剛形成的輻射激波(20 ns 之前)上游存在一個溫度低于1 eV 的預加熱區(qū), 這便是波前結構, 波前隨著激波溫度的降低很快消失.隨著激波進一步向外擴張, 激波陣面逐漸與高溫區(qū)表面分離, 形成了溫度的雙臺階狀分布(1 μs以后的溫度分布形狀), 高溫區(qū)和激波的溫度也不斷下降.

根據氙氣對外的輻射能流可計算其等效溫度Teff, 選取光學厚邊界(將345 nm—420 nm 波段的輻射吸收系數從最外層向內累加到1 的球層)為氣體發(fā)光位置, 在這一位置將不同波段的對外輻射譜能流按譜積分得到總能流S, 利用可以計算Teff, 其中σ是斯提芬-玻爾茲曼常數.圖4 給出了等效溫度的時間變化曲線, 30 ns—1 μs 的數值抖動是高溫區(qū)邊界的空間網格分辨率不夠導致的.等效溫度的變化過程表現出兩個明顯的峰值和一個極小值, 其中第一個峰值的等效溫度(約3 eV)要高于第二個峰值(約1.5 eV).下面詳細研究激波演化過程中氙氣的輻射特征及對應的物理過程.

圖5 給出了高溫氙氣對外輻射光譜的時間演化過程, 計算的波長覆蓋100—3600 nm.從圖5 看出不同波段輻射能流的演化規(guī)律也與等效溫度類似, 即存在兩個輻射強度的峰值, 并在2 μs 左右所有波段的輻射都減弱.為了研究氙氣輻射光譜的特征, 圖6 截取1 μs (圖6(a))和6 μs (圖6(b))兩個時刻氙氣的對外輻射光譜(實線), 并與相同等效溫度下的黑體輻射光譜(虛線)作比較.對比發(fā)現這兩個時刻高溫氙氣的輻射光譜明顯偏離黑體輻射光譜, 且不同時刻輻射光譜的形狀也不一致.由于計算模型使用多群近似, 每個能群內的輻射強度被平均化處理, 計算的輻射光譜不能顯示氣體輻射的精細結構, 但相對黑體輻射光譜的偏離是明顯的.

圖4 氙氣中輻射激波等效溫度( T eff , eV)的時間演化曲線Fig.4.Time evolution of the effective temperature ( T eff , in eV) of radiative shock in xenon.

圖5 輻射激波不同波段輻射能流的時間演化, 圖像亮度代表輻射能流大小Fig.5.Time evolution of radiation flux at different wavelength of the radiative shock, and the radiation flux are represented by the brightness of the figure.

為了研究氙氣輻射光譜的形成機制, 圖7 給出了不同時刻中心高溫區(qū)表面(虛線)和激波表面(點線)的黑體輻射光譜, 并與激波的對外輻射光譜(實線)作比較, 圖8 則給出了這些時刻380—460 nm波段輻射吸收系數的空間分布.在輻射激波擴張過程中, 隨著高溫氙氣不同位置溫度和密度的變化,這些位置的光學厚度也發(fā)生改變.氙氣的對外輻射光譜受輻射激波不同位置的輻射和吸收效應共同影響而偏離黑體輻射光譜, 并表現出觀察到的強度變化.

圖6 不同時刻氙氣的對外輻射光譜與相同等效溫度下的黑體輻射光譜的比較 (a) 1 μs; (b) 6 μsFig.6.Spectrums radiated from radiative shock in xenon and spectrums of blackbody radiation at the same effective temperature at different time: (a) 1 μs; (b) 6 μs.

剛形成的激波(10 ns)位于高溫區(qū)的外表面,激波表面和內部的溫度很高, 是光學厚的, 只有激波表面的輻射可以透出.此時激波的波前對輻射有較強的吸收(圖8(a)), 使對外輻射強度弱于激波表面輻射(圖7(a)).一段時間后(20 ns), 波前吸收變得很弱, 激波表面的輻射可以透出, 對外輻射光譜與激波表面的黑體輻射光譜接近(圖7(b)).在40 ns 時刻, 激波對400 nm—1500 nm 波段的吸收弱于其他波段, 使得這段輻射略強于激波表面輻射, 對外輻射能譜也偏離黑體輻射光譜(圖7(c)).這段時間由于波前吸收的減弱, 激波對外輻射呈增強趨勢, 對應圖4 等效溫度第一個峰值的上升段.注意在不同階段, 氙氣對于高能光子(波長在100 nm以下)的吸收始終很強.

激波在之后的擴張過程中逐漸變得透明, 使得激波后的部分輻射也可以透出, 對外輻射開始強于激波表面的輻射(圖7(d)和圖7(e)).到2 μs時刻,激波除對高能光子外幾乎變得完全透明, 高溫區(qū)表面的輻射可以透出(圖8(f)), 使得對外輻射能譜與高溫區(qū)表面的黑體輻射光譜接近(高能光子除外),氣體發(fā)光面過渡到高溫區(qū)表面(圖7(f)).這段時間激波區(qū)域的輻射強度隨溫度的下降而不斷降低, 對外輻射強度整體呈減弱趨勢, 對應圖4 等效溫度第一個峰值的下降段.

圖7 高溫區(qū)表面(虛線)、激波表面(點線)的黑體輻射光譜與氙氣對外輻射光譜(實線)的比較.(a)?(i)代表不同時刻的計算結果, 輻射能流統(tǒng)一投影到激波表面以避免球幾何發(fā)散帶來的影響Fig.7.Comparison among the blackbody spectrums at the surface of high-temperature-core (dash line), at the surface of the shock(dotted line), and the radiation spectrum outward-emitted by the radiative shock (solid line).(a)?(i) represent the result at different time, and all the radiation fluxes are projected on the shock surfaces to eliminate the spherical spreading effect.

圖8 380?460 nm 波段的輻射吸收系數隨位置的分布, 虛線為由外向內累積光學厚度為1 的位置.(a)?(i)代表不同時刻, 時間與圖7 相同Fig.8.Variation of absorption coefficient at 380?460 nm with position, the dash line mark the position where the optical depth equal to 1 from outside-in.Panel (a)?(i) are the same time as in Fig.7.

隨著溫度的降低, 高溫區(qū)表面也開始變得透明(圖8(g)), 內部的部分輻射可以透出, 使對外輻射強于高溫區(qū)表面的輻射(圖7(g)).隨著高溫區(qū)內部的輻射不斷透出, 氙氣的對外輻射又開始增強(圖7(h)), 對應圖4 等效溫度第二個峰值的上升段.最后, 高溫區(qū)溫度降低導致的輻射減弱效果超過輻射吸收系數降低導致的輻射增強效應, 對外輻射又開始降低(圖7(i)), 對應圖4 第二個峰值的下降段.

為了對比不同壓強下氙氣的對外輻射特征, 本工作還調整參數計算了輻射激波的光學演化過程.研究發(fā)現, 壓強在0.45—8 atm范圍內, 激波發(fā)光特性演化的主要過程是類似的, 只是特征時間與發(fā)光強度有所差異.所以, 以上的研究結果具有普遍意義.

4 總 結

本文使用Zinn 模型研究了氙氣中X 射線加熱產生輻射激波的發(fā)光特性及其演化過程.首先構建了一套氙氣的輻射不透明度與狀態(tài)方程參數, 并通過對比實驗數據, 驗證了模型和氙氣參數的合理性.

模擬研究發(fā)現, X 射線加熱的氙氣中形成了熱波、輻射激波、以及激波與內部高溫區(qū)的分離等現象.伴隨著這些過程, 氙氣的對外輻射表現出兩個亮度峰值和一個亮度極小值.且不同時刻的輻射光譜分布經常偏離黑體輻射光譜.這是由氣體中不同位置的輻射吸收系數隨時間的演化造成的.

剛形成的輻射激波位于高溫區(qū)的外表面, 激波陣面以內是光學厚的.此時氙氣的發(fā)光面位于激波表面, 但波前對輻射有強烈的吸收.隨著波前吸收的減弱, 氣體對外輻射逐漸增強并達到第一個峰值.之后, 隨著激波的擴張, 其溫度不斷降低, 對外輻射也不斷減弱.溫度降低也導致激波逐漸變得透明, 波陣面后的部分輻射可以透出.激波變得完全透明時, 發(fā)光面過渡到內部高溫區(qū)表面, 高溫區(qū)內部仍然是光學厚的.隨著溫度的進一步降低, 高溫區(qū)表面也逐漸變得透明, 內部的輻射可以透出, 氣體輻射強度隨之增強并到達第二個峰值.最后, 溫度降低導致內部輻射強度衰減, 氣體對外輻射強度再次降低.以上便是氙氣中輻射激波的演化過程,在這些過程中氙氣對波長100 nm 以下的輻射始終是不透明的, 而對400 nm—1500 nm 波段的吸收經常弱于其他波段.

本文的研究結果表明, 輻射激波的對外發(fā)射光譜并不嚴格符合黑體輻射光譜, 在黑體輻射假設下的計算中需要考慮這一偏離造成的誤差.在實驗診斷中, 由于高溫氣體的發(fā)光位置是動態(tài)變化的, 根據外部光學測量反推的物理量只是氣體不同區(qū)域對外輻射疊加后的等效值.

猜你喜歡
實驗模型
一半模型
記一次有趣的實驗
微型實驗里看“燃燒”
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
做個怪怪長實驗
3D打印中的模型分割與打包
NO與NO2相互轉化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 9丨情侣偷在线精品国产| 成人国产免费| 青草国产在线视频| 成人午夜免费观看| 97在线免费视频| 成人va亚洲va欧美天堂| 99热这里只有精品免费| 亚洲中文精品久久久久久不卡| 九九热精品在线视频| 幺女国产一级毛片| 欧美有码在线| 福利一区在线| 国产特级毛片| 青青草国产在线视频| 欧美日韩国产在线人成app| 国产成人精品第一区二区| 亚洲视频无码| 国产毛片高清一级国语 | 国产福利一区二区在线观看| 日韩精品免费在线视频| 欧美日韩一区二区三区四区在线观看| 国产午夜福利片在线观看| 国产午夜福利在线小视频| 欧洲一区二区三区无码| 成色7777精品在线| 精品伊人久久久久7777人| 四虎亚洲精品| 国产午夜无码专区喷水| 成年人免费国产视频| 免费无码又爽又黄又刺激网站 | 亚洲高清无码精品| 第一页亚洲| 一级全黄毛片| 亚洲日本在线免费观看| 白浆免费视频国产精品视频| 人妻精品全国免费视频| 91在线播放国产| 国产成人1024精品下载| 成人综合在线观看| 日韩a级毛片| 欧美一级视频免费| 国产精品福利社| 亚洲 成人国产| 国产波多野结衣中文在线播放| 91免费国产在线观看尤物| 欧美色香蕉| 国产极品美女在线观看| 国产福利免费在线观看| 免费xxxxx在线观看网站| 免费无码又爽又刺激高| 日本精品视频| 亚州AV秘 一区二区三区| 久久情精品国产品免费| 亚洲天堂啪啪| 亚洲无码高清视频在线观看| 91精品国产麻豆国产自产在线| 第一区免费在线观看| 久久毛片免费基地| 国产午夜福利片在线观看 | a在线亚洲男人的天堂试看| 中文字幕伦视频| 久久精品国产精品国产一区| 国产成人高清精品免费| 无码国内精品人妻少妇蜜桃视频| 亚洲日韩国产精品综合在线观看| 婷婷综合色| 亚洲AV无码乱码在线观看代蜜桃| 日韩A∨精品日韩精品无码| 欧美综合成人| 毛片网站观看| 亚洲国产无码有码| 欧美一级特黄aaaaaa在线看片| 伊人久久福利中文字幕| 99久久无色码中文字幕| 国产精品无码一区二区桃花视频| 在线不卡免费视频| 黄色一级视频欧美| 狠狠色狠狠综合久久| 成人va亚洲va欧美天堂| 亚洲av无码专区久久蜜芽| 欧洲亚洲欧美国产日本高清| 成人va亚洲va欧美天堂|