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

用于激光尾波加速的彎曲毛細管內氣流運動的模擬研究*

2023-10-06 07:04:34趙月琪崔佩霖李建龍李博原祝昕哲陳民劉振宇
物理學報 2023年18期

趙月琪 崔佩霖 李建龍 李博原 祝昕哲 陳民 劉振宇?

1) (上海交通大學機械與動力工程學院,上海 200240)

2) (上海交通大學物理與天文學院,上海 200240)

高壓放電充氣毛細管可產生等離子體通道,用于激光尾波加速.為探究尾波級聯加速所使用毛細管內的氣體流動及分布規律,本文建立了基于標準k-ε 模型的彎曲毛細管內氣體流動計算模型.以氦氣為工質,對彎曲毛細管內可壓氣體流動過程進行數值模擬,分析了不同結構、充氣背壓、充氣口位置對毛細管內氣體密度分布及速度場的影響.結果表明: 雙側對沖彎曲毛細管在充氣口之間管段具有較為穩定的氣體密度分布,充氣口附近氣體密度波動隨充氣口與毛細管兩端距離的增大而減小;在“直+彎”結構的級聯加速毛細管中,負責電子注入的直通道口徑會對彎管內氣體密度分布造成影響,當電子注入通道口徑小于150 μm 時,彎曲毛細管內氣體流動受到直通道的影響較小,可作為級聯結構中的電子束導引通道.

1 引言

超短超強激光在低密度等離子體中可激發等離子體尾波,該尾波具有超高的加速梯度(高出傳統射頻電子加速器3 個數量級),使得激光尾波加速有望成為新一代緊湊型電子加速的方案,成為近年來國內外研究的熱點[1-4].激光尾波加速需要克服聚焦光束長距離傳輸存在的自然散焦問題.利用高壓放電的毛細管可產生具有橫向拋物型密度分布的等離子體通道,可導引激光,大大延長聚焦激光的傳輸距離,從而增大激光尾波場的加速長度,提高電子束能量.近年來,等離子體通道導引的高能激光尾波加速被廣泛關注[5-8].Gonsalves 等[8]通過在直徑800 μm、長度20 cm 的充氣毛細管中放電電離氫氣,優化充氣背壓及放電延時后成功獲得7.8 GeV 的電子束.想要進一步提升電子加速能量,則需要利用多束激光驅動的尾波級聯加速.為此,Steinke 等[9]通過等離子鏡反射引入新的激光,利用等離子體透鏡聚焦上級出射電子束到下一級尾波場,初步實現兩級激光級聯的尾波電子加速.然而,由于等離子體透鏡的色散性質,電子束在兩級之間的耦合效率被大大降低(耦合率僅達到3.5%).2018 年,Luo 等[10]提出一種利用“直+彎”的毛細管結構產生等離子體通道,實現電子級聯加速的方案;其優點是避免使用具有色散性質的等離子體透鏡,有望提高電子束在兩級之間的耦合效率.對該類特殊毛細管內氣流的分析是設計該類等離子體通道的前提,目前尚未見報道.

在毛細管導引的激光等離子體尾波加速方案中,為確保激光能夠在等離子體中穩定長距離傳輸,獲得沿毛細管軸向均勻分布的等離子體密度尤為重要.而毛細管內等離子體由其內部氣體通過高壓放電產生[11],因此,充氣過程中毛細管內氣體流動狀態與最終等離子體密度分布直接相關.其難點之一是設計出氣體工質沿軸向均勻分布的毛細管結構,而復雜形狀管內的微尺度氣體流動特性直接影響其密度分布,需要進行詳細的數值模擬研究.目前,針對微尺度氣體流動數值模擬的方法主要有計算流體動力學方法[12-16]、格子-玻爾茲曼方法(lattice Boltzmann method,LBM)[17-20]和直接模擬蒙特卡羅(direction simulation Monte Carlo,DSMC)方法[21-23]等,其中關于彎曲邊界的研究已有較多文獻報道.Yan 等[15]分別用二階速度滑移邊界條件和分形幾何學方法描述氣體的邊界滑移并構建隨機粗糙表面,采用計算流體動力學方法研究了微通道內氣體在速度滑移和隨機表面粗糙度耦合作用下的流動特性,發現隨機表面粗糙度對微尺度下氣體的邊界滑移有顯著影響.閆晨帥和徐進良[16]采用 SSTk-ω低雷諾數湍流模型研究了加熱條件下超臨界壓力CO2在水平圓管內三維穩態流動和傳熱特性,通過實驗數據對比驗證了數值模型的可靠性和準確性,分析了熱流密度和質量流速對水平圓管內超臨界壓力CO2流動換熱的影響.Dai 等[18]提出了一種改進的曲線邊界處理方法,采用二階速度滑移條件的多松弛時間LBM 模擬微氣體流動,較好考慮了晶格節點與物理邊界之間實際偏移的影響,但未考慮氣體可壓的影響.Shariati 等[23]使用DSMC 算法改變微通道內隨機產生的圓形障礙物曲折度研究了多孔微通道中的氣體輸運行為,對比不同氣體類型流體動力學特征發現氣體本身性質會影響介質的表觀滲透性.LBM和DSMC 方法可較準確描述微尺度氣體流動特性,但存在模擬區域局限、模擬結構單一、計算耗費資源大的缺陷,不適合整段彎曲毛細管內氣體流動過程的全流程模擬.

本文使用商業ANSYS FLUENT 17.0 模擬軟件進行模擬,借助標準k-ε湍流計算模型,對不同結構及工作條件下毛細管內的氦氣工質流動過程進行數值模擬.探究不同充氣方式(單側直沖、雙側對稱)、充氣背壓(68950–172375 Pa)、充氣位置(距毛細管兩端3–12 mm)及不同結構(彎曲毛細管、“直+彎”毛細管)對毛細管內氣體流動的影響,揭示影響管內氣體密度波動規律,并獲得毛細管內氣體流動過程中的速度場、壓力場及管內中心軸線上氣體密度的變化規律.所得研究結果有望為基于曲率漸變彎曲毛細管的激光尾波級聯加速實驗提供理論指導及技術支持.

2 模型建立

2.1 幾何模型

如圖1 所示,放電毛細管裝置由兩端充氣口、彎曲毛細管和兩端銅電極組成,用于放電產生等離子體的氦氣由兩端充氣口注入,流向毛細管內以及兩側真空腔,并由裝置兩端電極連通高壓放電產生等離子體.由于毛細管為對稱性較好的圓管,綜合考慮計算耗時和計算精確性,此處將毛細管簡化為二維計算域,并在毛細管出口兩端增設真空計算區域以更為準確地模擬工質氣體由入口到真空腔的流動特性.其中毛細管長30 mm,內管徑0.5 mm.為獲得高壓充氣下較為理想的氣體流動過程,充氣段毛細管管徑設置為0.3 mm,過大的充氣段毛細管管徑導致充氣口與毛細管形成較為狹窄的流動通道,高速氣體沖擊至毛細管壁面致使更多氣體流入真空腔,中間管內流速增大,充氣口位置處形成較大渦旋.毛細管兩端設置真空計算區域長12 mm,寬10 mm.

圖1 放電毛細管結構示意圖Fig.1.Schematic diagram of discharge capillary structure.

2.2 數值模型

彎曲毛細管簡化為二維流體計算區域,其控制方程如下所示:

連續性方程為

能量守恒方程為

式中,ρ為流體密度,ui和uj為流體速度,p為流體壓力,μ為流體動力黏度,μt為湍流黏度,kf為導熱系數,cp為流體定壓比熱容,Prt=0.85 為湍流普朗特數,Tf為流體溫度.

對于穩態條件下彎曲毛細管內氣體流動過程,管道充氣口至出口部分為湍流狀態,靠近毛細管出口處管內流速最高可達923 m/s,根據25 ℃下氦氣動力黏度1.98×10-5Pa·s、管內氣體密度0.15 kg/m3計算得到雷諾數Re=ρVd/μ=3496,其中,ρ為流體密度,V為流體速度,d為圓管直徑,μ為流體動力黏度.標準k-ε模型比較適合用于求解高雷諾數湍流的情況,其輸運方程為

式中,K為湍流動能,ε為湍流動能耗散率,為速度分量,PrK和Prε為普朗克數,SK和Sε分別為K方程和ε方程的源項,Cε1,Cε2,Cε3為經驗系數.其中,

式中,Cμ為湍流常數,δij為克羅內克爾記號,為應變張量,βT為熱膨脹系數,Tˉ 為平均溫度,σT為湍流普朗克數,γ為定壓比熱與定容比熱之比,R為氣體常數.

對于具有較大流速(達到甚至超過當地音速)或者較大壓力變化的流體,其密度在流動過程中發生明顯變化,進而對速度場、溫度場及壓力場產生一定程度的影響,此時需要將流體作為可壓縮流體處理.通過馬赫數可以判斷是否需要考慮可壓縮效應:

根據理想氣體狀態方程,連續性方程、動量方程和能量方程分別可表示為

式中,T為氣體溫度;cV為定容比熱容;Rρ為ρ密度下氣體常數;ρ為可壓縮流體的密度,可由下式表示為

其中,pop+p為絕對壓力,Mw為氣體摩爾質量.

3 模擬條件及求解方法

3.1 網格劃分及參數設置

圖2 為彎曲毛細管結構的網格劃分,由于計算區域涉及彎管,因此采用結構化與非結構化耦合網格處理,有利于計算域數據傳遞的穩定性和收斂性,綜合考慮模擬精度與計算資源,選取管內氣體流速作為參考參數進行網格無關性驗證.結果如表1 所列,網格數達70628 后,管內氣體流速與前計算值的相對誤差減小到5%以內,最終選取70628 個網格作為毛細管流動模擬的網格數.毛細管內氦氣工質采用理想氣體模型,物性參數設置如表2 所列.

表1 網格無關性驗證Table 1. Grid independence verification.

表2 物性參數表Table 2. Physical parameters.

圖2 彎曲毛細管結構網格劃分及邊界設置Fig.2.Mesh division and boundary setting of curved capillary structure.

3.2 邊界條件及初始條件

模擬操作壓力設置為一個大氣壓(101325 Pa).

1)入口邊界: 壓力入口(pressure inlet),入口壓力103425 Pa,入口溫度為室溫(300 K).

2)出口邊界: 壓力出口(pressure outlet),考慮到出口為真空腔區域,真空泵壓力示數為5.5×10-4Pa(絕對壓力),因此模擬中出口壓力設置為-101324.99Pa,出口溫度為室溫(300K).

3)毛細管邊界: 毛細管壁面速度采用無滑移壁面,設置為絕熱壁面.

4)初始條件: 設置氣體計算域與真空腔的初始溫度為300 K,初始速度為0 m/s.

3.3 求解方法

放電毛細管內氣體流動數值模擬針對68950–172375 Pa 氣體注入壓力,距邊緣3–12 mm 充氣口位置,單/雙側對沖及“直+彎”毛細管結構進行研究.采用更適合可壓縮氣體基于密度的求解器對毛細管內氣體流動過程進行穩態模擬;使用標準k-ε模型,近壁面處理選用標準壁面方程,采用隱式格式進行計算,在精度滿足條件下進行合適的參數設置,以提高穩定性及收斂性.采用基于單元體的最小二乘法進行網格梯度離散,流動差分為二階迎風格式,湍流動能和湍流耗散率為一階迎風格式,庫朗數設置為5 以調節計算的穩定性和收斂性,湍流動能和湍流耗散率均設置為0.8,湍流黏度設置為1.連續性方程及氣體流速收斂精度為10-4,動量方程收斂精度為10-6,以此作為每次迭代計算收斂的判據.

4 結果與討論

4.1 模型可靠性驗證

取實驗充氣背壓即初始毛細管初始充氣壓力為4000 Pa、管內氣體流速為280–350 m/s 條件下測得的毛細管中心軸線上等離子體密度與毛細管中心軸線上氣體數密度進行數值模擬驗證,如圖3 和圖4 所示,分別為實驗和模擬密度趨勢對比圖與管內氣體流速分布云圖.其中,實驗等離子體密度分布由測量的等離子體發射特殊譜線的空間分布和碰撞展寬反演得到[7].

圖3 實驗等離子體密度與模擬氣體密度對比圖 (a)實驗等離子體密度;(b)模擬毛細管中心軸線上氣體密度Fig.3.Comparison of experimental plasma density and simulated gas density: (a) Experimental plasma density;(b) simulated gas density on the central axis of the capillary.

圖4 模擬計算域氣體流速分布Fig.4.Gas flow velocity distribution in the simulated computational domain.

等離子體由毛細管兩端電極對管內氣體高壓放電產生,管內氣體密度作為等離子體產生基礎條件可一定程度反映高壓放電后管內等離子體密度波動趨勢.由圖3 可知,模擬毛細管中心軸線上氣體數密度與實驗測得毛細管中心軸線上等離子體密度具有相同的變化趨勢,在充氣口附近密度呈明顯下降的趨勢,表現出激光導引段等離子體密度(氣體密度)不穩定,流入毛細管的氣體在兩個充氣口之間形成密度穩定段,與實驗測得的激光導引段管內等離子體密度分布相近[24].由圖4 可知,氣體在兩個充氣口之間的流速穩定在200 m/s.因此本數值模擬具有較好的可靠性,能夠較準確地預測管內氣體流動分布情況,為進一步在實驗上獲得穩定的激光導引提供數值預測.

4.2 充氣方式及充氣壓力對毛細管中心軸線上氣體密度影響

圖5 模擬了單/雙側充氣,壓力分別為68950,103425,137900,172375 Pa 時毛細管內中軸線上的氣體密度曲線圖.相較于雙側充氣,單側直沖毛細管內氣體密度波動較為劇烈,兩端充氣口位置出現較大的密度波動,中軸線上氣體密度并不完全遵循隨充氣壓力升高單調升高的規律,雙側對稱充氣毛細管在充氣位置中間段氣體密度分布較為穩定,形成平穩流動區域.圖6 給出了當充氣背壓為68950 Pa 和137900 Pa 時,在單/雙側充氣結構條件下充氣口位置附近的流線圖,氣體從充氣口注入毛細管通道時,由于充氣口與毛細管形成突擴結構,導致流動膨脹,流動邊界層在交界點處發生分離,由于分離后毛細管左右兩側壓力不同,導致自由剪切層發展過程中向出口處偏移,形成圓弧狀,左側剪切層與毛細管壁面底部發生碰撞,右側剪切層與毛細管壁面頂部發生碰撞,形成新的附面層,并在碰撞前的區域形成回流區,回流區形成的氣體渦旋導致充氣口位置的氣體密度降低.對比圖6(a)和圖6(b)與圖6(c)和圖6(d)不同壓力、不同充氣方式下流線圖發現,單側直沖充氣方式下,增大充氣壓力導致充氣口處氣體流速大幅增大,對彎曲毛細管壁形成沖擊后氣流一分為二,分別流入中間管段和管外真空腔,由于流入中間管段氣體流速較大,氣體對毛細管壁面形成較為明顯的沖擊,導致了中間管段內氣體流動不穩定,管內氣體流動愈加不穩定,而雙側對沖充氣方式下增大充氣壓力后,雙側進氣氣流在毛細管處形成交匯大部分氣體流入真空腔,中間管段保持著相對較低的管內流速,管內氣體流動較為穩定,可以觀察到此時毛細管充氣位置處已形成較大渦旋,進一步增大充氣壓力會加劇管內氣體流動的不穩定性.具體表現為: 在單側直沖結構毛細管充氣口位置附近形成的氣流渦旋較大,使得該區域氣體密度產生明顯波動;而雙側對沖結構毛細管充氣口位置附近的渦旋較小,氣體密度波動也較弱;然而隨著充氣背壓的增大,單/雙側充氣毛細管內氣體流動都愈加不穩定,毛細管內的氣體密度均出現明顯波動.模擬結果與實驗測得的不同充氣方式及充氣壓力下激光導引段管內等離子體密度分布趨勢相近[24].

圖5 不同充氣方式和充氣壓力下模擬管內氣體密度(a)單側直沖;(b)雙側對沖Fig.5.Gas density in the simulated tube under different inflation methods and inflation pressures: (a) One-side inflation;(b) double-side inflation.

圖6 不同充氣壓力下,不同充氣方式管內氣體流線圖 (a) 68950 Pa 單側直沖;(b) 68950 Pa 雙側對沖;(c) 137900 Pa 單側直沖;(d) 137900 Pa 雙側對沖Fig.6.Gas streamlines in the simulated pipe with different inflation methods: (a) One-side inflation under 68950 Pa;(b) double-side inflation under 68950 Pa;(c) one-side inflation under 137900 Pa;(d) double-side inflation under 137900 Pa.

4.3 充氣位置對毛細管中心軸線上氣體密度影響

圖7 為103425 Pa 充氣背壓下,毛細管中軸線上氣體密度分布圖(充氣口位置距毛細管兩端距離分別為3,6,9,12 mm).由圖7 可知,隨充氣口位置與毛細管兩端距離增大,管內氣體穩定區域長度不斷減小,在充氣口位置附近出現較明顯的氣體密度波動,波動幅度隨充氣口位置距毛細管兩端距離的增大而減小.隨著充氣位置距毛細管兩端距離增大,氣體從充氣管流入毛細管的壓降損失減小,毛細管內平均流速也隨之減小,位置距離毛細管兩端越近的充氣口位置氣壓損失更為嚴重,相較于充氣口距毛細管兩端12 mm 的氣體流動,氣體更多流向管外真空腔,導致中間管內密度相對較小,充氣口中間段氣體密度呈現平臺型密度分布.圖8 為當充氣口與毛細管兩端距離為12 mm 時,位于左右兩側充氣口之間的氣體流線圖,由于氣流運動路徑不同、與毛細管接觸部位的結構不同,導致了不同的壓降,表現為氣體從彎曲側流向直段.圖9 進一步給出了在相同距離下,雙側對稱充氣毛細管內的壓力云圖.由計算結果可知,氣體在進氣管道上壓強損失非常小,所以毛細管內氣壓值近似等于充氣氣壓,在靠近毛細管出口處,因為管內外壓強差導致毛細管內氦氣加速向外擴散,因此密度很快下降到與管外相同.因此在進氣口之間產生較均勻的氣體密度,通過調整充氣位置,可以實現對等離子體通道長度的調控.

圖7 不同充氣口位置毛細管中心軸線上氣體密度Fig.7.Gas density on the central axis of the capillary at different gas filling positions.

圖8 充氣位置距毛細管兩端12 mm 時,充氣口之間流線圖Fig.8.Streamline diagram between the left and right inflation port at 12 mm from both ends of the capillary.

圖9 充氣位置距毛細管兩端12 mm 時,充氣口之間壓力云圖Fig.9.Pressure cloud diagram between the left and right inflation port at 12 mm from both ends of the capillary.

4.4 電子注入通道口徑對彎曲毛細管內氣體密度分布的影響

為實現基于曲率漸變彎曲毛細管的激光等離子體尾波級聯加速,在上述彎曲毛細管的基礎上增加了實現上級電子注入的直通道[10].在不造成明顯激光能量損失的前提下,需要穩定且高效地將第二束激光脈沖導引到第二加速級中軸線上,而電子注入通道需要在彎管段與激光導引通道合并.第二加速級的結構如圖10 所示.

圖10 “直+彎”放電毛細管結構示意圖Fig.10.Schematic diagram of the “straight+curved” discharge capillary.

圖11 和圖12 為不同電子注入通道口徑條件下,彎曲毛細管中軸線上氣體密度曲線以及充氣口位置處的流線分布.由結果可得,當電子注入通道口徑較小時,氣體流動過程與雙側對稱充氣彎曲毛細管內的氣體流動過程類似,彎管中軸線上氣體密度擾動較小,充氣口位置處由于流體膨脹形成的渦旋導致該區域氣體密度波動較明顯.這一現象在100 μm 和150 μm 的電子通道口徑條件下表現得不太明顯,當電子注入通道口徑大于200 μm 時,如圖11、圖12(d)和圖12(e)所示,由于電子注入通道口徑的增大,彎管內氣體流動受到顯著的影響,充氣口處直管與彎管連接結構處氣體流動連通進一步加大,充氣口位置處亦形成較大渦旋,增大了該處管內中心軸線上氣體密度的擾動,直管彎管內氣體流動的連通也導致此處氣體密度下降,最終導致毛細管內氣體分布不均.因此,考慮管內氣體流動狀態與氣體密度均勻性,100 μm 和150 μm的電子注入口徑較適合“直+彎”結構實現電子在激光等離子體尾波中的級聯加速.

圖11 不同電子注入通道口徑下彎曲毛細管內的氣體密度分布 (a)直管徑50 μm;(b)直管徑100 μm;(c)直管徑150 μm;(d)直管徑200 μm;(e)直管徑300 μmFig.11.Gas density distribution in curved capillary tubes with different electron injection channel diameters: (a) 50 μm diameter;(b) 100 μm diameter;(c) 150 μm diameter;(d) 200 μm diameter;(e) 300 μm diameter.

圖12 不同電子注入通道口徑下彎曲毛細管內的氣體流速分布圖 (a)直管徑50 μm;(b)直管徑100 μm;(c)直管徑150 μm;(d)直管徑200 μm;(e)直管徑300 μmFig.12.Gas velocity distribution in curved capillary tubes with different electron injection channel diameters: (a) 50 μm diameter;(b) 100 μm diameter;(c) 150 μm diameter;(d) 200 μm diameter;(e) 300 μm diameter.

5 結論

本文針對上海交通大學新型激光尾波加速級聯方案,基于標準k-ε模型建立彎曲毛細管內氣體流動計算模型,對氦氣工質在曲率漸變彎曲毛細管內的流動過程進行模擬分析,獲得了彎曲毛細管內氣體密度分布情況,與實驗測量等離子體電子密度分布結果具有相同的變化趨勢,能夠較準確地預測管內氣體流動分布情況.基于該模型分別對單側直沖、雙側對沖和“直+彎”級聯加速結構毛細管內氣體流動過程進行模擬,得到以下結論:

1)相較于單側直沖彎曲毛細管結構,雙側對沖彎曲毛細管在左右兩側充氣口之間的氣體密度波動較小,氣體流動更為平穩,可產生較為穩定的等離子體密度通道.

2)在雙側對沖彎曲毛細管中,相同充氣背壓下毛細管兩個進氣口之間形成了較均勻的氣體密度分布;進一步研究結果表明,可通過控制充氣口的位置獲得不同長度的較均勻等離子體密度分布.

3)“直+彎”級聯加速毛細管結構中,電子注入通道口徑對彎管內氣體密度分布造成影響,電子注入通道口徑較小時毛細管內絕對壓力較低,充氣口與彎管之間的壓差較大,將導致彎管內氣體流速較高,使彎管中氣體密度波動變大;進一步研究表明,100 μm 和150 μm 的電子注入通道口徑較適合應用于“直+彎”級聯加速毛細管結構設計.

本文所得研究成果有望為基于曲率漸變彎曲毛細管的激光尾波級聯加速實驗提供理論指導及技術支持[25,26].

主站蜘蛛池模板: 国产第一色| 国产精品观看视频免费完整版| 特级aaaaaaaaa毛片免费视频 | 色爽网免费视频| 成人一级免费视频| AⅤ色综合久久天堂AV色综合| 精品国产黑色丝袜高跟鞋| 亚洲综合精品香蕉久久网| 亚洲精品免费网站| 亚洲六月丁香六月婷婷蜜芽| 亚洲成aⅴ人片在线影院八| 久久综合亚洲鲁鲁九月天| 亚洲最黄视频| 日韩福利在线观看| 国产激情无码一区二区三区免费| 黄色网页在线观看| 精品国产免费第一区二区三区日韩| 最新国产精品鲁鲁免费视频| 精品国产免费观看一区| 国产性生大片免费观看性欧美| 五月婷婷综合色| 中文字幕乱妇无码AV在线| 六月婷婷精品视频在线观看| 天天摸夜夜操| 玖玖精品在线| 免费在线观看av| 国产熟女一级毛片| 国产成人精品男人的天堂下载| 国产欧美精品专区一区二区| 欧美日本二区| 日本草草视频在线观看| 亚洲国产成人在线| 鲁鲁鲁爽爽爽在线视频观看| 婷婷色一二三区波多野衣 | 欧美成人影院亚洲综合图| 国产午夜精品鲁丝片| 2024av在线无码中文最新| 欧美一区精品| 国产办公室秘书无码精品| 孕妇高潮太爽了在线观看免费| 国产日本欧美在线观看| 精品三级在线| 亚洲成网站| 国产网站一区二区三区| 54pao国产成人免费视频| 小13箩利洗澡无码视频免费网站| 国产丰满大乳无码免费播放| 99热国产在线精品99| 91在线国内在线播放老师| 五月激情综合网| 国产激情无码一区二区APP| 日韩美毛片| 亚洲精品日产AⅤ| 一本色道久久88| 茄子视频毛片免费观看| 亚洲一区二区日韩欧美gif| 韩国福利一区| 国产精品护士| 国产福利在线观看精品| 黄色一级视频欧美| 国产日韩欧美成人| 国产精品吹潮在线观看中文| 欧美精品亚洲精品日韩专区va| 国产区福利小视频在线观看尤物| 久久久久夜色精品波多野结衣| 99爱在线| 老司机午夜精品网站在线观看 | 亚洲精品国偷自产在线91正片| 成人午夜亚洲影视在线观看| 五月婷婷综合网| 激情综合网激情综合| 欧美日本中文| 噜噜噜久久| 激情综合网激情综合| 亚洲无码高清一区二区| 亚洲欧美日韩另类| 青草视频久久| 污污网站在线观看| 久青草网站| 日韩精品亚洲一区中文字幕| 国产极品粉嫩小泬免费看| 免费一级成人毛片|