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

乏燃料剪切機氣固兩相流分析及結構優化

2020-10-24 01:41:24王時龍董建鵬周科源
原子能科學技術 2020年10期
關鍵詞:模型

黃 輝,王時龍,*,馬 馳,周 杰,董建鵬,周科源

(1.重慶大學 機械傳動國家重點實驗室,重慶 400044;2.中國原子能科學研究院 放射化學研究所,北京 102413)

乏燃料是在核反應堆中使用完畢后卸出的核燃料,國際上有一次通過長期地質貯存和閉式循環回收再利用兩種處理策略,我國確定實行后者,回收乏燃料中大量有用的核資源、分離長壽命放射性核素,結合嬗變技術,減少最終地質處置量,提高核能整體經濟性,從而保證核能的可持續發展,同時減少環境輻射污染。閉式循環回收再利用的一個重要環節就是采用以錒系元素萃取分離為主工藝的核燃料后處理,該流程就是對核燃料進行破解和溶解,進而進行萃取分離。燃料組件的破解是第一步,技術路線較多,目前普遍采用機械法剪切設備將核燃料處理成符合溶解需求的碎料。剪切機作為該設備的核心機械,設計時需考慮多方面的需求,如剪切時會產生大量強放射性核燃料芯塊粉塵,設計時需確保放射性粉塵盡可能保留在剪切機設備腔體內,防止擴散出污染熱室,導致無法進入熱室對其中的設備進行維修,必須將這些粉塵通過密封、氣流走向控制等設計首端進行嚴密隔離,并確保大部分乏燃料芯塊從剪切機流道內排入化學溶解器,保證放射性包容的同時,提高后處理的回收率。通過研究剪切機內污染性粉塵顆粒的分布和運動規律來降低剪切機內粉塵沉積水平是新型剪切機設計的重點。

傳統剪切機的設計并未涉及粉塵顆粒對乏燃料剪切機設計的影響以及氣固流邊界條件對粉塵顆粒運動分布的影響[1-4],因而通過傳統的剪切機設計方法[5-7]不能設計出滿足降塵溶塵要求的剪切機。剪切機內粉塵顆粒分布運動問題屬于流體動力學范疇,但一些研究流體動力學的方法,尚未在剪切機的設計上得到應用,且這些流體動力學方法在應用時,往往需要根據具體對象的特性和待研究的物理量等的不同來進行選擇[8-11],一般不能直接相互借鑒經驗。本文擬對乏燃料剪切機內氣固流的分布運動規律進行數值模擬研究,包括數值模擬方法的選取、流道結構和邊界條件對氣固流分布運動的影響,進一步對剪切機的結構和邊界條件進行優化,最終達到降低剪切機內粉塵沉積水平的目的。

1 模型建立和網格劃分

1.1 乏燃料剪切機實體模型建立

本文乏燃料剪切機的設計相比于傳統剪切機的設計,增加了通排氣和降塵溶塵的設計要求,其所剪切的乏燃料是由上、下兩段組成的材料。根據剪切機的功能要求以及設計理論,建立了剪切機幾何實體模型,如圖1所示。氣流運動的路徑大致為:外部通入的氣流分別從1個乏燃料進氣口、2個空氣噴嘴進氣口、1個箱體后區側壁進氣口進入到剪切機內,然后經過剪切處的氣流通道進入分料器溜槽,最后從溜槽排出后進入溶解器[12]。溶解器主體呈扁平圓柱形,其內部有溶解輪和溶解液,剪切機溜槽與溶解器的中心相連接,剪碎的乏燃料進入溶解器后,溶解輪圍繞溶解器中心步進式旋轉,使得碎料分散到溶解液中。由此剪切機內的污染性粉塵顆粒在氣體的攜帶作用下進入溶解器中被溶解吸收,達到降塵溶塵的目的。

1.2 剪切機內流場實體模型建立

由于剪切工位在剪切乏燃料時的模型最復雜,故本文選取的建模對象是剪切機剪切乏燃料時形成的相對結構,對其建立流體域模型,且考慮到本文研究的重點是剪切機內的流體域,在確保關注部位分析精度的前提下,對剪切機進行如下假設和簡化:1) 忽略幾何模型中小的圓角及倒角;2) 相對固定的結構視為一個整體;3) 忽略剪切機箱體外部的機構附件;4) 忽略剪切機內部機構的內部附件;5) 剪切機內產生的粉塵顆粒用Al2O3球形顆粒代替。簡化后的幾何模型如圖2所示。在ANSYS DM中對簡化幾何模型的內流場進行抽離,得到剪切機內流場實體模型,如圖3所示。

圖1 剪切機幾何實體模型Fig.1 Solid model of shearing machine

圖2 剪切機簡化幾何模型Fig.2 Simplified solid model of shearing machine

圖3 剪切機內流場模型Fig.3 Flow field model of shearing machine

1.3 內流場實體模型網格劃分

采用ICEM CFD軟件對流體域模型進行網格劃分,選用非結構化網格。同時采用Patch Conforming Method對較重要區域的網格進行細化和加密,以提高數值計算的精度。剪切機內流場的網格劃分如圖4所示。

圖4 剪切機內流場的網格劃分Fig.4 Meshing of flow field

1.4 內流場數值模型建立

1) 氣相數值模型建立

(1) 流動狀態確定。本文研究對象的環境溫度為25 ℃、大氣壓力為101 kPa。考慮到剪切阱附近的結構尺寸較大且大部分氣流處在剪切阱附近,選取過剪切阱大截面處的氣流來評估剪切機內氣流的流動狀態。本文的氣流管道為矩形截面,選取雷諾數(Re)的計算公式為:

Re=4vR/μ

(1)

R=ab/2(a+b)

(2)

其中:v為氣流平均速度;R為矩形通流截面的水力半徑;a和b分別為矩形流通截面的長度和寬度;μ為氣流的動力黏度。

根據剪切機實際工作要求,在截面處測得v=3.2 m/s;矩形截面的長度為1.25 m、寬度為0.7 m,R=0.224 4;查表得氣體的動力黏度μ=1.844 8×10-5Pa·s。計算得Re=155 698,大于2 300,可知剪切機內的氣流處于湍流狀態,所以采用湍流模型進行模擬計算。

(2) 湍流模型選取。從計算能力和精度看,Standardk-ε以及RNGk-ε模型是目前常用的2種模型。Standardk-ε模型中,湍動能方程k是精確推導計算,而湍流耗散率輸送方程ε建立在經驗公式上。RNGk-ε模型建立在嚴格的統計技術上,與Standardk-ε模型相比,其處理高、低雷諾數的復雜湍流的能力更強[11],本文選用RNGk-ε模型。

(3) 數值模擬方法選取。湍流數值計算有3種常用的數值方法:直接數值模擬法、大渦模擬法和雷諾平均法。雷諾平均法解決了直接數值模擬法計算量大的問題,避免了瞬時N-S方程的復雜非線性造成的精確描述三維流場以及時間特性相關全部細節的困難,其通過時均化的N-S方程將瞬態脈動量表現出來,對于解決工程實際問題,湍流的整體效果(即平均流場的變化)才是最重要的,而了解每個瞬態流場的全部細節對于解決工程實際問題幫助不大[13]。因此,本文采用雷諾平均法建立氣相數值模型。

(4) 氣相控制方程。由于剪切機內氣體流速遠低于聲速,氣相的動力黏度μf以及密度ρf可認為不變,忽略熱交換和體積力影響,得到氣相控制方程。

連續性方程為:

(3)

動量方程為:

(4)

其中:τ為氣相應力張量;β為相間動量交換系數;ε為耗散率;t為時間;uf為顆粒所在位置的流場速度;up為顆粒的運動速度。

2) 顆粒相數值模型建立

(1) 兩相流中顆粒的運動及受力分析。由牛頓第二定律可知,兩相流中任一顆粒的受力滿足mpdup/dt=∑F,mp為顆粒質量,∑F為顆粒所受外力的合力,包括顆粒受到的重力、浮力、斯托克斯阻力、壓力梯度力、附加質量力、貝塞特力、馬格努斯力、薩夫門力、曳力、泳力等。通過分析和計算,有些力數值很小,對數值模擬幾乎無影響,可忽略不計,本文主要考慮顆粒受到的重力和斯托克斯阻力的作用。

(2) 顆粒與壁面的相互作用。顆粒與壁面碰撞時,相互作用形式主要有反彈、粘附和逃逸。不同作用形式分別對應不同的壁面函數。在本文研究中,可近似認為顆粒與壁面之間的碰撞結果為反彈,且假設壁面無滑移,故考慮采用Scalable壁面函數[14]。

(3) 顆粒相基本運動方程。本文以湍流模型為基礎,采用在拉格朗日坐標系下進行顆粒運動追蹤的隨機軌道模型,用瞬時運動方程和隨機方法模擬湍流。由于模擬的粒徑和雷諾數相對較小,所以用經典模型來描述。為追蹤粒子,引入以下假設:追蹤顆粒是具有相同直徑和密度的球體;顆粒密度遠大于流體的密度,忽略浮力、壓力梯度力及其他較小的力;雖然剪切瞬間剪切位置粉塵顆粒產生量較密集,但考慮到顆粒會快速擴散及本文研究對象是剪切機的整個氣固兩相流而非剪切位置處的局部流體,此時顆粒相的體積濃度很小,可視為稀相,忽略粒子間的相互碰撞作用;忽略流體中的顆粒體積(非常小),顆粒的存在不影響氣體壓力;忽略溫度變化對流場的耦合影響及空氣含水量對空氣參數的影響。

基于上述假設,只考慮顆粒的重力和斯托克斯阻力。根據牛頓第二定律,單個粒子的運動方程可表示為:

|uf-up|(uf-up)/8+mpg(ρp-ρf)

(5)

在式(5)中代入顆粒阻力系數CD及顆粒相對雷諾數Rep的公式,顆粒的運動方程可改寫為:

dup/dt=3μCDRep(uf-up)/

(6)

up=uf+aτp/f+(up-uf-aτp/f)·

exp(-fΔt/τp)

(7)

xp=xf+(μf+aτp/f)Δt+

τp(μp-μf-aτp/f)(1-exp(-fΔt/τp))

(8)

(4) 氣固兩相耦合模型的選取。對于顆粒相為稀相的氣固兩相流,顆粒相和流體相之間的相互作用可分為單向耦合和雙向耦合兩種。本文顆粒體積濃度較低,一般情況下采用單向耦合即可滿足研究要求[15]。

2 數值求解方法和邊界條件

2.1 對流方程的離散

CFX在計算流體域模型前,要先對計算區域進行離散化處理。目前,有3種應用最廣泛的離散方法:有限元法、有限體積法和有限差分法[16],其中有限體積法綜合了其他2種方法的優點[17],求解方便,且求解精度滿足本文研究對象的要求。因此,選用有限體積法進行模擬計算。

CFX在進行非穩態求解時,有兩種精度選擇:一階向后歐拉以及二階向后歐拉。一階格式一般用于流體流動方向與網格成一條線的情況,此時求解的數值耗散相對較低;當流體流動方向與網格不在一條線時,對于各種多邊形網格,一般均采用二階格式,以此降低離散的誤差,提高求解精度,特別是當流體流動很復雜時,二階格式能較一階格式獲得更好的計算結果[18]。所以,本文采用二階離散格式來對流場模型求解。

2.2 流場迭代求解方法

本文采用的CFX可對各種流動進行質量守恒和動量守恒的求解。氣固兩相流耦合的問題主要是氣相和固相之間速度和壓力的耦合,由此需要解決以下2個問題:1) 在使用常規網格以及中心差分方法進行離散時,離散后的動量方程不能判斷壓力場是否合理;2) 在動量守恒方程中,壓力一階導是以源項表示的。針對第1個問題,可通過采用交錯網格來解決,而對于第2個問題,現階段最可靠的解決方案是SIMPLE算法。該算法的核心是“猜測-修正”原則,通過不斷的迭代,最終得到滿足要求的收斂解。該算法可求解交錯網格上的壓力場,進而對動量方程求解,是目前計算不可壓流場問題的主要方法。所以本文采用SIMPLE算法來計算。

2.3 邊界和輸入條件

根據剪切機實際工作要求,測量并通過收斂得到剪切機的進出口邊界條件:進氣口1氣流流速為1.000 m/s,進氣口2為1.000 m/s;進氣口3為3.740 m/s;出氣口的相對壓力為-1 100 Pa。輸入條件包括流場中氣體和顆粒的物性參數。流場中氣體的物性參數如下:材料為空氣;溫度為25 ℃;密度為1.169 kg/m3;動力黏度為18.448 μPa·s。根據剪切機實際工作情況,剪切阱附近是產生粉塵顆粒的主要位置,故在氣固兩相流模型中選取乏燃料剪切阱附近的3個位置(1、2、3)進行顆粒注入。顆粒的物性參數如下:材料為Al2O3、球狀、密度為3 700 kg/m3、直徑為50 μm、初始速度為1.000 m/s、數量為300×3。

3 數值計算結果及分析

3.1 氣相模擬計算結果

本文模型為流固耦合,采用二階差分和SIMPLE算法,求解量較大,綜合CFX與Fluent的優劣對比,采用CFX更適合本文的計算。本文使用CFX 19.0進行計算。

考慮到剪切機內流場的形態復雜,無法采用規則的平面或區域來描述流場中氣流的形態,本文選取3個平面(經過剪切阱,剪切阱區域如圖5所示)來描述流場中氣流的流動,3個平面分別為:平面1(z=3.245 m)、平面2(y=0.08 m)、平面3(x=-0.452 66 m)。

模擬得到的剪切機內平面1、2、3上氣流速度和湍動能的分布如圖6、7所示。在平面1上,速度最大的區域在剪切刀附近,為12.78 m/s,湍動能也主要分布在剪切刀附近,分布范圍為0.224 8~0.899 0 m2/s2。在平面2上,速度最大的區域為剪切箱體后區靠近偏心滾筒的位置,分布范圍為4.126~16.500 m/s,湍動能也主要分布在該區域(0.374 0~1.494 0 m2/s2)。在平面3上,速度最大的地方出現在剪切阱附近,為19.04 m/s,且可見大面積的渦流出現,湍動能為0.426 3~1.705 0 m2/s2,主要分布在進料口中程處以及溜槽處。

圖5 剪切阱及3個參考平面Fig.5 Shear well and three reference planes

剪切機內氣流壓力、速度分布流線圖如圖8所示。從圖8可看出,剪切機內氣流最大壓力出現在剪切口上方,最大速度出現在刀具附近。流線圖顯示氣流在剪切室后區形成大面積的渦流,阻礙了氣流的運動,使得粉塵顆粒隨氣流一直在剪切箱體內打轉而無法順利排出剪切箱體。分析可知,氣流旋渦的形成主要是由進氣口3的氣流運動造成的。綜上所述,對剪切阱附近的結構以及剪切箱體后區進氣口位置進行優化非常迫切。

圖6 平面1、2、3上氣相速度矢量圖Fig.6 Vector diagram of gas velocity on plane 1, 2 and 3

圖7 平面1、2、3上的湍動能分布Fig.7 Turbulent kinetic energy distribution on plane 1, 2 and 3

圖8 剪切機內氣流壓力和氣流速度流線圖Fig.8 Pressure streamline and velocity streamline of gas

3.2 結構改進及模擬

由前述模擬結果和分析可知,箱體結構可優化位置有兩處:剪切阱附近的結構、箱體后區側壁進氣口位置。由于本項目缺乏剪切阱附近結構的設計要求,無法對其作出改動,故本文選擇對箱體后區側壁進氣口位置進行優化。優化方案為:將剪切機箱體后區側壁進氣口去掉,改為從箱體下方的乏燃料上段溜槽口進氣。為實現上段溜槽口進氣均勻,可在上段溜槽口下方接一漏斗形進氣嘴。改進后得到的內流場模型及進出氣口分布如圖9所示。

改進進氣口3位置后,進行數值模擬,得到平面1、2、3上氣流的速度和湍動能分布,如圖10、11所示。分析圖10、11可知:改進進氣口3位置后,在平面1上,氣流速度最大處出現在溜槽的中部,約為23.45 m/s,湍動能主要分布在溜槽兩側壁面處,約為0.177 1~0.711 6 m2/s2,氣流速度最大的位置及湍動能的分布和大小較之前有所改善,可降低對剪切阱附近結構的破壞;在平面2上,剪切阱處的氣流流速較大,約為3.576~14.300 m/s,湍動能主要分布在剪切阱附近,約為0.184 4~0.737 5 m2/s2;在平面3上,乏燃料剪切處周圍的氣體流速較大,約為3.960~15.84 m/s,湍動能主要分布在溜槽前后壁面,約為0.275~1.100 m2/s2,與改進前的剪切處相比,可在一定程度上起到降低刀具損耗速度的作用。比較改進前后3個平面上的氣流速度分布以及湍動能分布可知,進氣口3位置改進后,同一平面上的速度和湍動能均有所下降,且其空間分布更有利于降低剪切刀的損耗,并在一定程度上延長剪切機的工作壽命。

圖9 改進后的內流場模型Fig.9 Improved flow field model

進氣口3位置改變后,通過CFX模擬,剪切機內氣流的壓力、速度分布的流線圖如圖12所示。分析圖12可知,改進進氣口位置后,氣流的形態有了很大改善,沒有了改進前的大面積旋渦,氣流路徑較為通暢,氣體從各進氣口進入后經過剪切室內部順利從出氣口排出,能有效攜帶剪切箱體內的粉塵顆粒進入溶解器。

圖10 改進后平面1、2、3上的氣相速度矢量圖Fig.10 Vector diagram of gas velocity on plane 1, 2 and 3 after optimization

圖11 改進后平面1、2、3上的湍動能分布Fig.11 Turbulent kinetic energy distribution on planes1, 2 and 3 after optimization

綜上,對進氣口3位置的優化是成功的,優化后剪切機內的氣流形態有了很大改善,能促進剪切機內粉塵顆粒的順利排出和降低剪切刀的損耗速度。

3.3 氣固兩相模擬結果

根據剪切機實際工作情況,剪切阱附近是產生粉塵顆粒的主要位置,故在氣固兩相流模型中選取乏燃料剪切阱附近的3處位置進行顆粒注入,注入位置如圖13所示。

剪切機內流場中3個平面上穩定后的顆粒體積分數分布的CFX模擬結果示于圖14,剪切機內氣流流線圖和顆粒運動軌跡示于圖15,顆粒運動動態模擬的開始-過程-結束3個時間點的截圖示于圖16。

分析圖14可知:平面1上的顆粒主要分布在剪切刀前方位置,平面2和平面3上顆粒主要分布在剪切阱附近。從圖15可看出,顆粒均流向溜槽一個出口,符合實際情況。但流體域內局部區域仍存在渦旋,不利于顆粒排出剪切機。從圖16可看出,顆粒注入流體域后隨氣流一起運動。在顆粒運動模擬開始時,每條流線上均有1個顆粒,這些顆粒從顆粒注射位置發出,沿各自的流線開始運動,此時剪切機內的顆粒數最多;在顆粒運動模擬過程中,各條流線上的顆粒沿各自的流線繼續運動,其到達的位置處于動態變化,此時一部分顆粒還在剪切機內運動,而另一部分已排出剪切機并進入溶解器,此過程為顆粒的擴散過程,剪切機內的顆粒數逐漸減少;在顆粒運動模擬結束時,各條流線上顆粒到達的位置已穩定,不再發生變化,此時顆粒的擴散結束,剪切機內的顆粒數最少。分析可知,大部分顆粒最終可順利排出剪切機,只有少部分顆粒留在剪切機內,如圖16中紅圈所標的顆粒即在剪切機內死角處形成沉積,不能排出剪切機。

圖12 改進后剪切機內氣流壓力和氣流速度流線圖Fig.12 Pressure streamline and velocity streamline of gas after optimization

圖13 剪切機內顆粒注入位置Fig.13 Particle injection location

圖14 平面1、2、3上顆粒體積分數分布 Fig.14 Distribution of particle volume fraction on planes 1, 2 and 3

圖15 剪切機內氣流流線圖(a)和顆粒運動軌跡(b)Fig.15 Airflow streamline diagram (a) and particles trajectory (b)

圖16 顆粒運動模擬截圖Fig.16 Screenshot of particle motion simulation

3.4 輸入條件優化

分析可知,剪切機內流場中顆粒相所占的比重Y與流道尺寸L、進出氣口的氣流流速v、流量Q和壓力p等有關,其函數關系可表示為Y=f(L,v,Q,p,O),其中,L、v、Q、p、O為獨立變量,O表示其他自變量。從純數學的角度看,這是多約束非線性規劃問題,可建立函數模型并通過遺傳算法進行求解。考慮到本文研究的剪切機對L、v、Q、p的限制,在其他條件限定不變時,允許進氣口1的流速v1可變,因此研究了v1對Y的影響。v1=16、17 m/s時,Y的值分別為1 777和6 094,異常大;v1=18 m/s時風級已達8級,為大風,對剪切機內部器件破壞性較大,故v1不再向后取值。剔除v1=16、17、18 m/s的3個點,通過擬合得到Y與v1的關系,如圖17所示。

圖17 Y與v1的關系Fig.17 Y vs. v1

分析可知,v1=9 m/s時,剪切機內顆粒物濃度水平達到最低值,故在實際工作中可將v1設定為9 m/s,以利于粉塵顆粒排出剪切機。

4 結論

1) 利用CFD軟件和數值模擬方法,研究了進氣口位置和輸入氣流流速對剪切機內固體顆粒運動及沉積的影響。進氣口位置改進后,剪切機內氣相的壓力和速度分布得到很大改善,大部分區域速度較均勻,氣流路徑較為通暢,避免了大渦流的出現;進氣口1的氣流流速為9 m/s時,剪切機內顆粒物的濃度水平達到最低值。利用同樣的方法,可研究氣體流量、壓力和流道尺寸等對剪切機內流場的影響。

2) 基于歐拉-拉格朗日方法的離散顆粒模型較適用于剪切機內稀相氣固兩相流的研究。剪切機內的氣流是典型的湍流,最大湍流動能出現在速度突變的位置,如剪切機內的死角或有渦流處。剪切機內顆粒物的運動與氣流形態密切相關,顆粒的分布與顆粒源的位置以及氣流形態密切相關,在顆粒源附近和有渦流的位置,顆粒濃度較高。

對剪切機內氣流形態和顆粒物分布運動的模擬研究對于優化剪切機結構和工藝有較強的指導作用,可為剪切機內氣流流動和粉塵控制提供理論上的支持。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 日韩精品亚洲人旧成在线| 国产剧情无码视频在线观看| 999福利激情视频| 为你提供最新久久精品久久综合| 超碰色了色| 嫩草在线视频| 欧美日韩在线国产| 亚洲视频影院| 免费看a级毛片| 天天干伊人| 91伊人国产| 99热这里只有成人精品国产| 99热线精品大全在线观看| 丁香婷婷综合激情| 欧美日韩国产系列在线观看| 国产精品毛片一区| 99视频只有精品| 国产在线高清一级毛片| 免费看美女自慰的网站| 毛片视频网| 国产欧美高清| 国产一区二区精品福利| 亚洲国产日韩在线成人蜜芽| 国产十八禁在线观看免费| yy6080理论大片一级久久| 99无码中文字幕视频| 波多野结衣一区二区三区四区| 日本AⅤ精品一区二区三区日| 9久久伊人精品综合| 欧美另类视频一区二区三区| 精品视频第一页| 久一在线视频| 国产欧美精品一区二区 | 国产成人高精品免费视频| 亚洲无码视频喷水| 国产玖玖玖精品视频| 尤物在线观看乱码| 97超碰精品成人国产| 综合五月天网| 国产在线自乱拍播放| 国产三级精品三级在线观看| 精品撒尿视频一区二区三区| 久久人妻xunleige无码| 欧美午夜久久| 国产成人狂喷潮在线观看2345| 欧美一区精品| 99热最新网址| 亚洲啪啪网| 日韩高清欧美| 日韩精品亚洲一区中文字幕| 欧美成人午夜影院| 国产人成在线视频| 国产菊爆视频在线观看| 97视频免费在线观看| 国产手机在线小视频免费观看| 国产亚洲精久久久久久久91| 狠狠色噜噜狠狠狠狠奇米777| 在线a视频免费观看| 日韩人妻无码制服丝袜视频| 99久久亚洲综合精品TS| 热re99久久精品国99热| 国产爽歪歪免费视频在线观看 | 国产视频a| 一区二区无码在线视频| 亚洲高清在线天堂精品| 国产精品xxx| 日本不卡在线播放| 亚洲天堂视频在线播放| 97久久免费视频| 国产SUV精品一区二区| 毛片在线播放a| 国产在线日本| 日韩在线网址| 九月婷婷亚洲综合在线| 欧美高清视频一区二区三区| 亚洲精品午夜天堂网页| 国产精品视频导航| 欧美乱妇高清无乱码免费| 91久久天天躁狠狠躁夜夜| 中文国产成人久久精品小说| 欧美亚洲另类在线观看| 老司国产精品视频91|