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

高重復頻率水冷Nd:YAG 激活鏡放大器的溫度特性*

2021-02-06 04:31:04肖凱博鄭建剛3蔣新穎蔣學君吳文龍嚴雄偉王振國鄭萬國3
物理學報 2021年3期
關鍵詞:模型

肖凱博 鄭建剛3) 蔣新穎? 蔣學君 吳文龍嚴雄偉 王振國 鄭萬國3)

1) (中國工程物理研究院激光聚變研究中心,綿陽 621900)

2) (中國工程物理研究院,高能激光科學與技術重點實驗室,綿陽 621900)

3) (上海交通大學,IFSA 協同創新中心,上海 200240)

為解決高重復頻率大能量激光放大器的熱管理問題,采用數值模擬與實驗分析的方法,對背面水冷Nd:YAG 激活鏡放大器的流體散熱進行了研究.基于低雷諾數k-ε 湍流模型,建立了流-固共軛傳熱多物理場耦合分析模型,對比分析了近壁面處理方法對流體流動、對流擴散和熱傳導過程及溫度分布的影響,分析研究了不同冷卻液流量和泵浦參數對流場特性、激光介質溫度和波前分布的影響.數值模擬表明: 激光介質的溫度分布與固液邊界層內的黏性作用密切相關,且冷卻液的熱擴散主要發生在100 μm 范圍內; 激光介質的熱沉積分布中心對稱,而溫度分布沿水流方向不對稱,最大溫升位于出水口端且基本保持不變; 增益介質前表面的溫度分布與介質的波前分布隨冷卻液流量非線性變化,而隨泵浦參數線性變化; 實驗結果與數值模擬符合較好.

1 引 言

二極管泵浦的高能重頻固體激光器因轉換效率高、亮度高、可靠性好和結構緊湊等優點而備受關注,同時其在工業加工、科學研究和國防軍事等領域也具有十分廣闊的應用前景.目前,輸出能量達到十焦耳甚至百焦耳、重復頻率運行的納秒激光系統在國內外得到了廣泛的研究[1?4].美國勞倫斯-利弗莫爾實驗室的Mercury 激光裝置采用疊片放大器結構,利用高速氦氣對介質表面進行主動冷卻,獲得了61 J/10 Hz 的輸出[5]; 采用相似的疊片結構和低溫氣冷Yb:YAG 片,英國STFC-RAL的DiPOLE100 激光裝置實現了最高的105 J/10 Hz 輸出[6]; 基于端面泵浦背面水冷的Yb:YAG激活鏡構型,法國LULI 實驗室的Lucia 裝置獲得了14 J/2 Hz 的激光輸出[7].這些激光裝置都是采用摻Yb3+的激光材料作為增益介質,盡管具有相對較低的熱負載率,但常溫下Yb3+離子是準三能級結構,閾值泵浦功率密度很高,對溫度波動異常敏感,飽和通量較高不利于能量提取,因而目前主要采用低溫冷卻提升材料的光譜和熱力學性質,但也增加了激光系統的復雜程度.相對而言,Nd:YAG 具有一些顯著的優勢,吸收截面和受激發射截面較大,便于提取能量,也能以合適的泵浦強度在室溫下運行,但對儲能和熱管理提出了一些新的要求.采用分布式背面水冷Nd:YAG 激活鏡構型,將增益和熱負載分散到多個放大器模塊中,有利于泵浦儲能和移除廢熱,近年來清華大學和本課題組均實現了12 J/10 Hz 的激光輸出[8,9].

然而,當前大能量脈沖激光器的重復頻率只能達到10 Hz 量級,難以滿足很多實際應用場景對更高重復頻率的需求[10,11],因此迫切需要研制百赫茲量級的重頻大能量激光器[12].在高泵浦功率密度下,重復頻率的提升將進一步加劇激光系統熱效應的累積,采用高效的冷卻結構、準確評估激光介質內的溫度分布是高重頻大能量激光器突破熱效應瓶頸,實現穩定可靠運行的重要途徑.在冷卻結構的設計和介質溫度評估中,大部分研究常采用對流換熱系數來簡化模型,而忽視了冷卻過程中冷卻液溫升和流場特性對介質溫度分布的影響[12?17].針對這一問題,李策等[18]通過解析表達,分析了流體直接冷卻板條激光器中對流換熱系數和冷卻液溫升與流速的關系,以及介質的溫度和應力分布.采用標準k-ε 或RNG k-ε 湍流模型,文獻[19?22]先后建立了熱-流-固多物理場耦合模型,對低溫氣冷疊片激光器、流體直接冷卻板條激光器、薄片激光器的流道結構和激光介質的溫度分布進行了研究.但上述湍流模型對黏性作用較強的固壁區采用壁面函數法作近似處理,并未細致評估壁面區的流動,且對流動分離過大或近壁面處于高壓下的流動也具有一定的局限性,因而會影響介質溫度分布的評估.低雷諾數k-ε 湍流模型能適應不同的雷諾數區域,對湍流區和近壁區均采用一套公式計算,具有較高的計算精度[23].目前,基于低雷諾數k-ε 湍流模型,對背面水冷Nd:YAG 激活鏡放大器進行熱管理的分析尚未報道,有必要進行深入研究,為高重頻水冷Nd:YAG 激活鏡放大器的優化設計和穩定運行提供指導.

圖1 水冷Nd:YAG 激活鏡放大器的模型 (a) 激活鏡放大器結構; (b) 流道結構Fig.1.Simulational model of water-cooled Nd:YAG active mirror amplifeir: (a) Configuration of the active mirror amplifier;(b) configuration of cooling channel.

本文主要研究了高重頻大能量背面水冷Nd:YAG 激活鏡放大器的溫度分布特性.基于低雷諾數k-ε 湍流模型,建立了激活鏡放大器流-固共軛傳熱多物理場耦合分析模型,利用COMSOL Multiphysics 有限元軟件對比分析了標準k-ε 和低雷諾數k-ε 湍流模型近壁面處理方法對流體流動、對流擴散和熱傳導過程及溫度分布的影響,分析研究了不同流量和泵浦參數下冷卻液的流場特性、激光介質溫度和波前分布,并開展了相關實驗研究,實驗結果與數值模擬符合較好.

2 物理與數學模型

2.1 水冷Nd:YAG 激活鏡放大器結構模型

二極管泵浦水冷Nd:YAG 激活鏡放大器的結構如圖1(a)所示,激光介質的背面沉浸在冷卻液去離子水中,側面和前表面由不銹鋼支撐固定; 高功率二極管陣列從Nd:YAG 背面泵浦,種子激光由其前表面注入經背面反射后雙程放大輸出.增益介質Nd:YAG 口徑為52 mm×32 mm,厚度為7 mm,摻雜濃度為1% (原子百分比),對泵浦光的吸收系數為3.0 cm–1,一面鍍1064 nm 高反膜和808 nm 增 透 膜,另 一 面 鍍1064 nm 增 透 膜 和808 nm 高反膜; 同時介質四周采用紫外固化膠粘合吸收系數為2.5 cm–1、寬為4.5 mm 的Cr4+:YAG包邊,以便抑制放大自發輻射和寄生振蕩.假定二極管陣列壓縮整形后輸出50 mm×32 mm 的平頂光束,最大泵浦強度為3 kW/cm2,并以90%的效率耦合進Nd:YAG 中被完全吸收.被吸收的泵浦光一部分直接轉化成熱,另一部分以反轉粒子數的形式存儲在激光上能級,然后又會以熒光或放大的熒光被包邊吸收而產熱[24].根據前期建立的泵浦和放大過程物理模型[25],假定Nd:YAG 的產熱率為37%,紫外固化膠的厚度為0.5 mm,考慮Cr4+:YAG 的吸收后,可得放大器激光介質內的熱源分布.為了簡化計算模型,忽略介質前表面的空氣對流冷卻及不銹鋼的傳導冷卻,其流道結構如圖1(b)所示,且以Nd:YAG 的背面中心為坐標原點; 介質背面與玻璃窗口之間的流道間隔為1.5 mm,冷卻液的初始溫度為20 ℃,其進入流道后先經過導流板,再流經Nd:YAG 背面進行冷卻.計算過程中的物性參數如表1 所列[14,26].

2.2 流-固共軛傳熱理論模型

在激光器泵浦和冷卻過程中,涉及流體流動、對流擴散和熱傳導等過程,其中流體的冷卻效果與流動狀態密切相關.一般認為湍流比層流的冷卻效果更好,其流動狀態可由雷諾數Re確定[23]:

式中,ρw,u和μ分別表示流體的密度、速度和動力黏度;Dh為水力直徑,對于矩形流道,為流道間隔的2 倍; 如果Re< 2300,流體可看作層流; 介于2300 與4000 之間,為過渡流;Re> 4000,可看作湍流.根據激活鏡放大器的流道結構和流速,本文主要針對湍流狀態進行分析.

用于冷卻的去離子水一般可視為不可壓縮黏性流體,其內部的流動狀態、對流擴散和熱傳導可通過連續方程、動量方程和能量方程來描述:

式中,t為時間;u為速度矢量,ui(i=1,2,3)分別是沿x,y和z方向的速度分量;p表示流體的壓強;xi(i= 1,2,3)分別表示x,y,z;kw和Cpw分別表示流體的熱導率和比熱容;T表示溫度;Qthw表示流體中的熱源.如果(4)式中的速度為0,則變為固體內的熱傳導方程:

式中,ρNd,kNd和CpNd分別表示固體的密度、熱導率和比熱容;QthNd表示固體中的熱源.

在流-固耦合邊界上熱耦合條件為

式中,nNd和nw分別為固體和流體耦合邊界的外法線矢量.在以往大多數研究中,通常采用對流換熱冷卻來簡化處理流-固耦合邊界,此時對流換熱系數h難以準確確定,從而影響計算的準確性.

式中,T∞為流體內溫度.

對于固壁區內的流動,由于分子黏性作用,流體速度將逐漸降低,此時Re數較低,湍流發展并不充分.然而,目前常用的標準k-ε 模型和RNG k-ε模型等都是針對充分發展的湍流才有效; 為了求解這一區域內的流動,這些模型引入了壁面函數法,采用一組半經驗的公式將壁面上的物理量與湍流核心區的相應物理量聯系起來,而不對黏性影響比較嚴重的區域進行有效地計算,因而并不能細致地評估壁面區的流動,從而也會影響流-固邊界層內的換熱[23].為此,本文采用低雷諾數k-ε 模型,使數值模型能適應不同的雷諾數區域,并結合流-固邊界層熱耦合條件,分析流體內的流動狀態及與固體之間的傳熱.

表1 計算過程中的物性參數Table 1.Physical parameters used in the simulation.

3 數值模擬與實驗驗證

3.1 溫度與流場分布

圖2 (a) 激光介質中橫向熱沉積分布; (b) 流道中心截面流速分布比較; (c) 標準k-ε 模型和(d) 低雷諾數k-ε 模型中激光介質及固液邊界層中的溫度分布Fig.2.(a) Transverse heat density distribution in the laser medium.(b) Comparison of flow velocity distribution across the center of the cooling channel.Temperature profiles in the laser medium and the solid and liquid boundary layers obtained by (c) the standard k-ε model and (d) low Re k-ε model.

根據圖1 所示的泵浦方式,在泵浦光傳輸方向上,泵浦光能量按吸收定律呈指數分布.假設放大器以50 Hz 的重復頻率運行,泵浦強度和脈寬分別為1.90 kW/cm2和240 μs,此時激光介質中橫向熱沉積分布如圖2(a)所示.從圖2(a)可知,激光介質橫向熱沉積分布呈中心對稱,Cr4+:YAG 吸收放大自發輻射光線后產生的熱量主要集中在邊界中心部分而在四個角上較少,且在邊界處的平均熱功率最大可達2.33×107W/m3,遠大于Nd:YAG中平均熱功率1.08×107W/m3.假設冷卻液在流道入口處初始溫度為20 ℃,入口總流量為10 L/min,出口參考靜壓為0,基于標準k-ε 模型和低雷諾數k-ε 模型計算所得的流道中心截面(x = 0)流速分布、激光介質及中心截面固液邊界層中的溫度分布,如圖2(b)—(d)所示.從圖2(b)可以看出: 由于導流板的水流勻化作用,冷卻液經過導流板后流速達到最大值,然后立即降速混合均勻,在流道末端,導流板也可減少湍流區域; 然而,當冷卻液進入激光介質背面轉接口時,流速變化較快,為了消除其對介質冷卻的影響,設計中延長了介質與兩端轉接口之間的流道距離以便介質背面水流更加均勻.對比標準k-ε 模型和低雷諾數k-ε 模型,兩種模型流道內的流速形貌基本一致,然而介質背面固液邊界層內標準k-ε 模型中流速相對更大,界面處甚至達到1.27 m/s,這主要是由于標準k-ε 模型忽視了黏性作用對邊界層的影響,壁面函數法也無法像低雷諾數k-ε 模型那樣細致地評估壁面區內的速度分布.由圖2(c)和圖2(d)可知: 在兩種模型中,激光介質內的最大溫升均位于介質表面包邊與Nd:YAG 邊界處,冷卻液流經激光介質吸熱后溫度會逐漸升高,從而導致介質內溫度分布與熱沉積分布不再一致,沿水流方向不再對稱; 然而由于邊界層內的黏性作用,在低雷諾數k-ε 模型中固液邊界層及介質內的溫升表現得更加明顯,標準k-ε 模型中由于邊界層內的流速較大從而導致縱向熱擴散很少,水流中的縱向溫升幾乎不存在.在低雷諾數k-ε 模型中,與Nd:YAG 相比包邊內的產熱更高且橫向熱傳導明顯,因而Nd:YAG 前表面的溫度分布沿水流方向將稍微降低、然后逐漸升高,其平均溫度達到49.62 ℃,溫差約為3.27 ℃; 而在后表面,由于水流換熱更加顯著,Nd:YAG 中的溫度將隨水溫升高而一直升高;當水流經激光介質時,固液邊界層內的溫升主要發生在100 μm 范圍內,其中界面的溫升主要來源于激光介質的傳熱,最大溫升可達10.85 ℃,此后熱擴散起主要作用,溫升降到1.2 ℃.

圖3 激光介質沿 (a) x 方向、(b) y 方向、(c) z 方向的應力分布; (d) 激光介質沿z 方向的形變分布Fig.3.Thermal stress distributions in the laser medium along (a) x direction,(b) y direction,(c) z direction; (d) z-deformation of the laser medium.

3.2 應力應變及波前分布

激光器運行時,激光介質內的溫度梯度將使介質在高溫區域產生的熱膨脹受到低溫區域的制約,形成不均勻的位移,從而產生熱應力.根據上述研究,將低雷諾數k-ε 模型獲得的激光介質溫度分布與固體力學進行耦合,可得激光介質的熱應力和形變分布,如圖3 所示.圖3(a)—(c)分別為激光介質沿x,y,z 方向的熱應力分布,可以看出,由于包邊內溫度梯度較大且前表面形變被約束,應力主要存在于包邊內,其中最大應力達到361 MPa,且沿x,y 方向的應力分布與溫度分布相似.從圖3(d)可以看出,由于紫外固化膠的膨脹系數較大且彈性模量較小,因而沿z 方向的形變略大于激光介質,最大達到3.37 μm; Nd:YAG 介質內沿z 方向的形變中心區域較為均勻,邊緣梯度較大,Nd:YAG介質內最大形變量位于后表面,達到為3.27 μm.

圖4 Nd:YAG 中熱致波前畸變(λ = 1064 nm) (a) 熱光效應引起的畸變; (b) 形變引起的畸變; (c) 總波前畸變Fig.4.Thermally induced wavefront distortions in the Nd:YAG slab: (a) Thermal-optical effect induced wavefront distortion;(b) the deformation induced wavefront distortion; (c) the overall wavefront distortion.

此外,介質內的溫度梯度也將導致光束通過放大器后光程差發生變化,從而引起波前畸變.對于激活鏡放大器,熱致波前畸變主要來源于四方面:1)激光介質由于熱膨脹產生的形變; 2)熱光效應導致介質折射率發生變化; 3)熱應力引起的彈光效應對折射率的影響; 4)激發態反轉粒子引起的電子作用對折射率的影響,其中后兩項的作用相對較小[13].根據上述所得的溫度、熱應力和形變分布可得光束通過Nd:YAG 介質后的熱致波前畸變,如圖4 所示.從圖4 可知,熱光效應引起的波前畸變與溫度分布基本一致,波前PV (peak-to-valley)值為0.78λ; 形變引起的畸變與熱光效應引起的畸變相反,對總波前具有較大的影響,且中心低、邊緣高,四周頂角熱畸變較大,約7.45λ; 由于熱光效應對形變的補償,Nd:YAG 中總波前畸變PV 值降到7.27λ,RMS (root mean square,RMS)值約為2.88λ.

3.3 流量對激光介質溫度和波前分布的影響

在散熱過程中,激光介質表面的對流換熱與流道的設計,流體的物性參數以及受熱條件等密切相關,其中流體流量是影響換熱的一個關鍵因素.在上述研究的基礎上,基于低雷諾數k-ε 模型,模擬分析了不同流量下激光介質的溫度和波前分布.假定冷卻液在流道入口處總流量從2.5 L/min 增加至22.5 L/min (間隔為2.5 L/min),放大器的運行條件保持不變,其計算結果如圖5 所示.從圖5 可知,隨著流速增加,Nd:YAG 前表面的平均溫度和溫差都將逐漸降低,但增益介質中最高溫度點始終保持不變,且當流量超過12.5 L/min 后,溫度變化逐漸平緩,溫度分布總體形貌基本一致; 同時,由于溫度梯度降低,波前畸變也相應減少,其PV 值由7.92λ 降低到7.13λ,最后趨于穩定.

3.4 泵浦參數對激光介質溫度和波前分布的影響

在放大器運行過程中,泵浦參數,如泵浦光強度和重復頻率等,將直接影響介質中的熱負載,從而影響放大器的熱效應.假定放大器以50 Hz 重頻運行,泵浦脈寬為240 μs,將泵浦電流從100 A 增加到390 A,相應的泵浦強度從0.63 kW/cm2遞增到2.95 kW/cm2,流量為10 L/min,則泵浦強度對Nd:YAG 溫度和波前的影響如圖6(a)和圖6(c)所示.隨著泵浦強度的增加,Nd:YAG 介質前表面平均溫度和溫差、波前PV 和RMS 值均線性增加,當泵浦強度達到2.95 kW/cm2時,平均溫度和溫差分別達到67.89 ℃和3.54 ℃,波前PV 和RMS值分別為12.17λ 和4.80λ.此外,假定放大器泵浦強度和脈寬分別為1.90 kW/cm2和240 μs,流量保持不變,改變放大器的運行頻率,則激光介質的溫度變化和波前分布如圖6(b)和圖6(d)所示.隨著運行頻率從10 Hz 增加到60 Hz,增益介質前表面的平均溫度和溫差、波前PV 和RMS 值也近似線性增加,當重復頻率低于35 Hz 時,Nd:YAG 前表面平均溫度和溫差分別在41 ℃和1.5 ℃以內,波前PV 和RMS 值分別在4.93λ 和1.95λ 以內.

圖5 流量對Nd:YAG 前表面溫度分布和介質波前的影響 (a) 中心截面(x = 0)處溫度分布曲線; (b) 平均溫度和沿水流溫差分布; (c) 波前分布Fig.5.Temperature field distributions and wavefront distortions of gain medium at different flow rates: (a) Temperature profiles across the center of gain medium; (b) mean temperature and temperature difference distributions; (c) wavefront distortions.

圖6 Nd:YAG 前表面平均溫度、溫差及介質波前隨(a),(c)泵浦功率和(b),(d)重復頻率的變化Fig.6.Average temperature and temperature difference on the front surface of Nd:YAG and the wavefront distortions of the gain medium at different (a),(c) pump intensities and (b),(d) repetition rates.

3.5 實驗驗證

為了驗證理論分析模型,采用紅外熱像儀對Nd:YAG 激活鏡放大器的表面溫度分布進行了測量,并與低雷諾數k-ε 模型的分析結果進行比較.放大器運行時,二極管泵浦電流為250 A,泵浦脈寬為240 μs,相應的泵浦強度為1.9 kW/cm2,通過觀察前表面的熒光分布來保證泵浦沉積的對稱性,同時為了避免種子激光對測量結果的影響,測量時放大器沒有通光.放大器在50 Hz 重復頻率下運行時,Nd:YAG 介質前表面實驗測得的溫度分布與數值模擬結果整體形貌基本符合; 但實驗測得的最高溫度點位于頂端右側,這主要是由于壓縮整形后泵浦源左右不對稱或二極管陣列巴條的泵浦強度不一致所致,如圖7 所示.此外,調節泵浦源的重復頻率,測量了放大器10—60 Hz 運行下Nd:YAG 前表面的溫度分布,其最高溫度點處沿水流方向的溫度分布曲線,前表面的平均溫度和溫差,如圖8 所示.從圖8 可知: 隨著熱負載的增加,Nd:YAG 前表面的溫度逐漸升高,且最高溫度點位置保持不變; 對于Nd:YAG 前表面的平均溫度,實驗測量和數值模擬的結果基本符合,且兩者之間的差異與頻率無關; 當頻率小于25 Hz 時,實驗測量與數值模擬的溫差基本一致,此后隨頻率的增加,兩者之間的差異逐漸變大,這主要是由于數值模型忽略了Nd:YAG 背面介質膜層的熱阻和不銹鋼支撐件的熱傳導所致.

圖7 50 Hz 重復頻率下Nd:YAG 前表面的溫度分布 (a) 數值模擬值; (b) 實驗測量值Fig.7.Temperature field distributions on the front surface of Nd:YAG operating at the repetition rate of 50 Hz: (a) Numerical simulation results; (b) experimental results.

圖8 (a) 實驗測量的Nd:YAG 最高溫度點處沿水流方向的溫度分布曲線; (b) Nd:YAG 前表面平均溫度和溫差的理論與實驗對比Fig.8.(a) Measured longitudinal temperature profiles through the peak temperature of Nd:YAG; (b) theoretical and experimental comparison of the average temperature and temperature difference on the front surface of Nd:YAG.

4 結 論

本文對高重頻大能量背面水冷Nd:YAG 激活鏡放大器的溫度分布特性進行了數值模擬與實驗研究.建立了背面水冷激活鏡放大器流-固共軛傳熱多物理場耦合分析模型,考慮到固液邊界層內的黏性作用,相比標準k-ε 湍流模型,低雷諾數k-ε 湍流模型對壁面的處理更加細致,可較精確地模擬固液邊界層內的流動狀態及介質的溫度分布.在1.9 kW/cm2和240 μs 的泵浦條件下,放大器以50 Hz 運行時,固液邊界層中冷卻液的熱擴散主要發生在100 μm 范圍內,其中界面的最大溫升可 達10.85 ℃,Nd:YAG 前 表 面 的 平 均 溫 度 為49.62 ℃,沿水流方向的最大溫差為3.27 ℃,波前PV 和RMS 值分別為7.27λ 和2.88λ.由于冷卻過程中,冷卻液吸熱后溫度逐漸升高從而導致介質溫度分布與熱沉積分布不再一致,沿水流方向也逐漸升高; 且Nd:YAG 前表面的平均溫度和溫差、介質波前PV 和RMS 值隨冷卻液流量非線性變化,而隨泵浦參數線性變化.實驗測量的Nd:YAG 前表面平均溫度與數值模擬結果基本符合,而隨著頻率增加,實驗測量與數值模擬的溫差具有一定的差異,這主要是由于數值模型忽略了Nd:YAG 背面介質膜層的熱阻和不銹鋼支撐件的熱傳導所致,所得結果對高重頻大能量水冷Nd:YAG 激活鏡放大器的優化設計和穩定運行具有一定的指導意義.對水冷Nd:YAG 激活鏡放大器的波前畸變進行實驗研究及補償將是下一步工作的重點.

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 九一九色国产| 国产精品林美惠子在线观看| 久久免费成人| 最新精品国偷自产在线| 国产人成乱码视频免费观看| 亚洲开心婷婷中文字幕| 日本伊人色综合网| 影音先锋丝袜制服| 午夜福利在线观看入口| 激情综合婷婷丁香五月尤物| 精品精品国产高清A毛片| 蜜桃视频一区二区| 欧美日韩第二页| 国产一级一级毛片永久| 日日拍夜夜操| 在线观看免费人成视频色快速| 亚洲综合激情另类专区| 人妻丰满熟妇AV无码区| 国产福利拍拍拍| 999福利激情视频| 国内毛片视频| 国产在线观看一区精品| 欧美成人午夜影院| 人人艹人人爽| 污污网站在线观看| 国产福利免费视频| 亚洲国产精品一区二区第一页免 | 久久香蕉国产线看精品| 一级毛片免费观看不卡视频| 久久久受www免费人成| 国产chinese男男gay视频网| 中文字幕亚洲精品2页| 国产精品丝袜视频| 免费又黄又爽又猛大片午夜| 久久综合丝袜长腿丝袜| 国产在线小视频| 欧美日韩另类国产| 欧美日韩国产综合视频在线观看 | 日韩在线观看网站| 爆乳熟妇一区二区三区| 日本www色视频| 国产91成人| 欧美日韩精品一区二区在线线| 亚洲成在线观看| 国产理论一区| 黄色网站不卡无码| 久久久久久国产精品mv| 国产毛片基地| 国产最新无码专区在线| 国产成人综合网在线观看| 色香蕉影院| 九九热在线视频| 无码不卡的中文字幕视频| 播五月综合| www.91中文字幕| 欧美在线综合视频| 伊人精品视频免费在线| 54pao国产成人免费视频 | 粉嫩国产白浆在线观看| 精品无码一区二区在线观看| 国产一级在线播放| 国产一区二区三区在线观看视频| 日韩国产一区二区三区无码| 成人福利在线视频| 大乳丰满人妻中文字幕日本| 亚洲精品成人片在线播放| 欧美日本在线观看| 夜夜操狠狠操| 91视频区| 91无码人妻精品一区二区蜜桃| 一区二区三区国产精品视频| 综1合AV在线播放| 国产成人狂喷潮在线观看2345| 日本欧美成人免费| 国产鲁鲁视频在线观看| 日韩av手机在线| 亚洲欧美成aⅴ人在线观看| 狠狠色噜噜狠狠狠狠色综合久| 在线视频97| 中文成人在线视频| 国内毛片视频| 国产精品自在自线免费观看|