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

氣體噴吹對富氧側吹爐內兩相流動影響的數值模擬

2022-03-29 03:13:48祝振宇周萍陳卓龍鵬張嶺
中南大學學報(自然科學版) 2022年2期
關鍵詞:風速區域

祝振宇,周萍,陳卓,龍鵬,張嶺

(1.中南大學能源科學與工程學院,湖南長沙,410083;2.長沙有色冶金設計研究院,湖南長沙,410001)

基于瓦鈕科夫側吹煉鉛技術,我國冶金企業聯合高等院校及研究院改進研發出富氧側吹煉鉛工藝[1]。自2010年以來,多家小型冶煉企業應用富氧側吹氧化熔煉-氧氣側吹還原熔煉工藝后,企業多項技術經濟指標都取得了顯著提高。

富氧側吹煉鉛工藝具有工藝流程與備料工序簡單、原料適應性強、爐床能力與一次出鉛率高、投資和生產成本低以及爐子壽命長等特點[2]。研究表明,富氧側吹煉鉛工藝的鉛回收率比富氧底吹煉鉛工藝的鉛回收率高1.38%[3]。近年來,富氧側吹熔煉工藝也逐步拓展應用至鉍、銻及工業廢料冶煉生產中[4-5]。

富氧側吹煉鉛工藝在我國的廣泛應用使得其相關的研究工作也得以迅速發展。2007年,在我國首屆熔池熔煉技術及裝備專題研討會上,賓萬達[6]針對我國富氧側吹煉鉛技術的大型工業化應用,提出了2個亟需解決的問題,即氧化區和還原區的連接問題以及鉛渣還原過程中如何實現“焦層過濾”的問題。這些問題成為后續研究的重點方向。張立等[7]對高鉛渣的還原過程進行了分析,認為盡管熔池內部熔體攪拌劇烈,但新生的鉛液滴在相互碰撞后聚合長大并沉降到噴槍以下的安靜區域,因此在鉛還原爐中不需要額外設置鉛渣分離區。李允斌[8]指出在富氧側吹爐氧化段爐體采用銅水套雖然造成了熔煉過程需要補熱的局面,但是這一措施對于噴槍和壁面耐火材料起到了良好的保護作用。白樺等[9-11]對氧氣側吹技術在鉛鋅共生氧化礦上的應用進行了研究,實現了冶煉廢渣的無害化處理。祁棟等[12]對富氧側吹爐處理廢舊電池鉛泥過程進行了初步的評述,認為在熔池內部產生劇烈攪拌的情況下,富氧側吹爐內的氣液固相間反應極快,爐床能力比其他爐型的高,適合用于處理廢舊電池鉛泥等。趙娜等[13]通過分析富氧側吹氧化爐的電除塵設備運行數據,認為氧化爐出口煙氣量、含塵濃度和硫酸濃度等波動較小,爐子整體運行穩定。經過多年的研究,富氧側吹煉鉛技術在廢舊材料和低品位精礦處理方面表現出巨大應用潛力,但富氧側吹煉鉛過程在整體能耗、泡沫渣[14]和煙塵率[1]處理方面仍具有進一步提高和完善的潛力,且有關其熔池內部熔體流動與傳熱傳質過程方面的系統性研究仍有待深入。

近年來,隨著計算機算力的飛速提升,數值模擬成為諸多領域中[15-16]不可或缺的重要研究手段,并在冶金多相體系的研究中得到了廣泛應用,包括對P-S 轉爐[17-19]、瓦鈕科夫爐[20-22]、煙化爐[23-24]和錳鐵脫碳爐[25-26]等冶金爐窯的生產過程解析,但針對富氧側吹爐的數值模擬研究較少。LIU等[27]對不同湍流模型在側吹熔池數值模擬中的適用性進行過綜合比較與分析,為進一步深入研究奠定了基礎。

為了進一步了解富氧側吹爐一次風噴吹方式和噴吹速度對熔池內熔體運動的影響,本文以鉛熔煉富氧側吹爐為對象,分別針對生產中常見的“對吹”和“錯吹”2 種噴吹方式,在風速175,225 和275 m/s 下開展數值模擬,以分析不同條件下富氧側吹爐內熔體的運動特點,從而為富氧側吹爐操作制度的優化提供重要參考。

1 數值模型的建立與驗證

1.1 物理模型

鉛熔煉富氧側吹爐的爐體結構大致可以分為上部的煙氣區域、中部的渣相區域以及底部的金屬相區域3個區域,爐體結構簡圖如圖1所示。爐體全長為9.75 m,高約6.10 m,渣相上部區域(熔池噴吹區域)的寬度為2.20 m。在渣相區域的爐體兩側分別布置有13對一次風口,呈高低交錯排布,其中低風口7 對,高風口6 對,各風口直徑均為32 mm。同側相鄰兩個風口之間相距0.75 m,高風口距離靜止渣液面的浸沒深度為0.500 m,低風口浸沒深度為0.625 m。在實際生產過程中,13對風口僅需要使用一半數量的風口即可滿足生產需要,因此,在生產操作中采用的主要噴吹方式有對吹和錯吹2種。

圖1 富氧側吹爐爐體結構示意圖Fig.1 Schematic of oxygen-enriched side-blown furnace

為提高計算效率,仿真計算中根據富氧側吹爐一次風口的分布規律選取部分爐體作為數值求解區域。鑒于爐體兩側的高低風口呈對稱交錯分布,建模時選取高、低風口各2 組構建幾何模型,如圖2(a)所示,其具體結構與尺寸如圖2(b)所示。在此部分爐體區域中,由于爐內渣液面以下區域的兩相流動過程為研究核心,且不關注頂部煙氣流動,因此,將煙氣擋板簡化省略。另外,采用混合網格對計算域分區進行網格劃分,其中熔池內部和上部氣相空間采用六面體網格,爐體中部使用四面體網格逐漸過渡(如圖2(c)所示)。鑒于風口附近氣體流速快、流動相對復雜的特點,對此區域內的網格進行局部加密處理。

圖2 數值模擬物理模型與網格劃分Fig.2 Geometry and mesh for simulation

1.2 數學模型

富氧側吹爐內兩相流動過程的數值模擬涉及流動、傳熱和傳質等過程,在對其進行多物理場耦合模擬時不僅需要考慮多相流模型和湍流模型的正確選擇,還需要考慮氣流的熱脹性等因素對兩相流動狀態的影響。富氧側吹爐內的流體介質包括氣相、渣相和金屬相3種,由于本研究的重點為一次風口處的渣相區域,故在數值模擬中將金屬相和渣相均視為熔體相而予以簡化。根據課題組前期對富氧側吹爐內熔體流動的數值模擬研究[27-28],選取Mixture 多相流模型和Realizablek-ε湍流模型對側吹爐內兩相流動過程開展數值模擬,相間曳力模型與滑移速度則分別采用常見的Schiller-Naumann 與Manninen 模型[29],相關數值模型的可行性驗證可參見文獻[27-28]。另外,本數值模擬采用非穩態的計算方式,在保證庫朗數合理的條件下,將計算的時間步長設置為0.1 ms。

數值模擬中求解的守恒方程組如下。

1)連續性方程:

式中,vm為混合相平均速度,m/s;ρm為混合相密度,kg/m3;t為時間,s。

2)動量方程:

式中:p為絕對壓力,Pa;μm為混合相的動力黏度,Pa·s;F為體積力,N;g為重力加速度,m/s2;vdr,k為任意相k的漂移速度,vdr,k=vk-vm。

3)能量方程:

式中:αk為第k相的體積分數;Ek為熱焓;keff為有效導熱系數;T為溫度;SE為源項;vk和ρk為第k相的速度和密度;n為相的個數。

1.3 邊界條件與物性參數

計算域內包括3類邊界:質量流量入口、壓力出口和壁面,分別對應圖2(a)中一次風入口深色面、計算域頂部紅色面和其余各淺色面,各邊界上的參數設置如表1所示。需要說明的是,由于采用結構化網格和非結構網格相混合、分區劃分網格的形式,且一次風口呈高低交錯排布,若按照慣常方法將計算區域前后端墻設置為周期性邊界,則會降低網格質量,因此,仿真計算區域中將此兩側端墻設置為壁面邊界,區域內雖然包括了4組風口,但僅以中間的2 組風口所在的區域(即風口所在中部的1.5 m爐長區域內)作為數據提取和分析區域,以減小壁面邊界條件對氣液兩相流動的影響。此外,質量流量入口處的湍流強度根據常溫下的一次風速及風口水力直徑計算得到,入口一次風氣流中氧氣與空氣的質量比為5∶1,但本數值模擬中暫不考慮氣液兩相間的化學反應;同時,由于壓力出口位于富氧側吹爐豎直煙道端部,根據生產中微負壓操作參數,將其設置為-10 Pa。

表1 數值求解邊界條件Table 1 Boundary conditions for simulation

仿真計算時所使用的氣相和熔體相物性參數如表2所示,熔體相的物性參數設置為高鉛渣的相關數值,其中高鉛渣的密度和比熱容按照其化學組成計算,高溫下高鉛渣的黏度采用文獻[30]中的數據。同時,考慮到氣液兩相間巨大的溫度差可能對氣體密度以及兩相流動過程的影響,將氣體設置為不可壓縮理想氣體,其密度僅隨溫度變化。在數值模擬的初始條件下,熔體相的溫度為1 473.15 K,一次風入口氣流的溫度為300 K。

表2 氣相與熔體相物性參數Table 2 Physical properties of gaseous and melt phases

1.4 模型驗證

為驗證網格尺寸和數學模型的合理性和準確性,采用如圖2所示的分區劃分網格的方式,在熔體區域分別采用3種不同尺寸的網格,獲取3套不同的網格體系以開展網格獨立性檢驗,如表3 所示。計算過程中,3套網格均采用相同的數學模型和邊界條件,硬件配置為4塊Intel Xeon Gold 6248處理器,共計80 核心,計算時長如表3 所示。為驗證模型的準確性,參考已有文獻[31-32]中的氣-液兩相系側吹實驗,并計算速度基本消失的無因次距離x′:

式中:x′為無因次距離;x為一次風口軸線上某一點至風口處距離,m;d0為一次風口的直徑,m。

氣-液兩相系側吹過程中,速度基本消失時的無因次距離x′為10~20[31-32]。由表3可知,網格1的無因次距離僅為8.98,而采用網格2 和網格3 的無因次距離與文獻[31-32]接近,且符合實驗所得結論。因此,綜合考慮計算時長等因素,后續計算均采用網格2尺寸開展計算,如圖2(c)所示。

表3 網格獨立性檢驗Table 3 Mesh independence study

2 數值模擬的工況設置與結果分析

富氧側吹爐內一次風的噴吹方式和風速是影響熔池內兩相流動過程的重要因素,也是爐體操作制度優化的重要方向。為了研究不同噴吹方式和噴吹速度條件下熔體的運動特點,在對吹和錯吹2 種噴吹方式下,進行風速分別為175,225 和275 m/s的數值模擬。

2.1 不同噴吹方式對熔體運動的影響

2.1.1 工況參數設置

一次風口作業方案示意圖如圖3 所示,其中,L1,L2,L3,L4,R1,R2,R3 和R4 分別為富氧側吹爐切片中的8個一次風口。錯吹與對吹工況開啟的風口位置及編號如紅色箭頭所示。錯吹與對吹工況參數設置如表4所示。

表4 不同噴吹工況下仿真參數設置Table 4 Simulation parameters in different gas injection modes

圖3 一次風口作業方案示意圖Fig.3 Schematic diagram of opening nozzles

2.1.2 數值模擬結果分析

由于仿真采用非穩態的計算方式,因此需要首先確定熔池內部兩相流動由開始噴吹到進入相對穩定狀態所需要的時間。為此,選取如圖4(a)所示熔體相靜止時所占灰色區域,計算熔體區域的平均速度:

式中:m為某空間區域內的網格節點數;vi為該區域內網格節點i處的速度,m·s-1;vˉ為該區域內的平均速度,m·s-1。

熔體區平均速度隨時間的變化曲線,如圖4(b)所示。從圖4(b)可見,在噴吹開始的0.75 s 內,熔體區域的平均速度變化較大,而在0.75~2.50 s 內,2 種噴吹方式下熔體區域的平均速度基本穩定在0.38 m/s 左右,波動幅度均小于±5%,可以認為在2 種噴吹方式下,熔池內部氣-液兩相運動狀態在此時段內都已基本趨于穩定。

圖4 熔體區域平均速度的變化曲線Fig.4 Changes of average velocity of melt zone

為深入了解熔池內不同高度位置處熔體相的速度和相體積分數等信息,取若干位置截面(如圖5(a)所示,其中,爐內熔體區域階梯處中部位置為坐標原點O,H為爐高方向,W為爐寬方向),計算不同高度截面上熔體相體積分數平均值,結果如圖5(b)所示。從圖5可以看出,對吹工況條件下的熔體相體積分數沿爐體高度方向上的變化規律與錯吹工況下的結果相似。但是在對吹工況條件下,渣液面高度以下的熔體相體積分數普遍低于錯吹工況的熔體相體積分數,而在渣液面高度以上位置處,其熔體相體積分數則相對高于錯吹工況的熔體相體積分數,這說明在對吹工況條件下,靜止液面以上區域內受一次風帶動而翻卷起來的熔體較多,由此推測當采用2組低風口對吹時,沖擊至熔體中的2股一次風氣流能夠使熔體翻卷的程度更高。

圖6 所示為t=2 s 時對吹和錯吹工況下不同爐體高度和寬度位置處熔體相的平均速度(面平均速度)的變化曲線。從圖6 可知,無論是沿爐高方向還是爐寬方向,對吹工況下的熔體相平均速度均高于錯吹工況下的結果。特別是從圖6(b)可以看到,當采用4 個低風口對吹時,靠近風口位置處(爐寬W=±0.75 m)的熔體相平均速度較為接近,均約為0.45 m/s,比錯吹工況中相同位置處的熔體平均速度分別約高0.15 m/s 和0.25 m/s。這說明采用低風口對吹時,在一定程度上有助于提高熔池內熔體的運動速度,結合圖5中對吹工況下熔體體積分數在靜止渣液面以上區域更大這一結果,說明對吹時熔體翻卷的程度更劇烈。而在圖6(a)中,在H=1.467 m位置,對吹工況下的熔體相平均速度有所下降,這是因為在該高度位置上,有部分熔體運動達到其最高點而處在即將回落的狀態,其熔體相速度較低。

圖5 熔池內不同高度位置及熔體相體積分數沿爐體高度的變化曲線Fig.5 Schematic of different heights and changes of melt volume fraction along height of furnace

圖6 t=2 s時熔體相平均速度沿爐高和爐寬方向的變化曲線Fig.6 Changes of melt average velocity along height and width of furnace at t=2 s

圖7所示為熔池區域的速度矢量圖,其中,圖7(a)中所示深色截面區域即為矢量圖所選截面區域,圖7(b)和(c)所示分別為錯吹工況和對吹工況下低風口該截面的速度矢量圖。從圖7可以看出,即使在錯吹工況下,一次風口上部位置熔體仍會出現類似于對吹工況下熔體對稱翻卷的狀態,說明左側位置熔體會受到相鄰位置處一次風的帶動,出現熔體向熔池中部翻卷的現象;但相較于對吹工況,錯吹工況下熔體的翻卷狀態受右側一次風影響更大,并略向左偏移。

圖7 熔體區域速度矢量圖Fig.7 Velocity vector in melt zone

為了更準確地掌握不同工況下熔體受一次風帶動而獲得的速度,分別計算如圖5(a)所示渣液面高度、低風口高度以及鉛液面高度以下空間區域內的熔體相平均速度,如表5所示。通過對比可以發現,在任意高度位置以下的空間區域內,對吹工況中的熔體平均速度均比錯吹工況時的高。這說明在相同的風量條件下,采用4個低位風口對吹時,氣流對熔體流動的影響更大,熔池中熔體的運動更為劇烈,因而在一定程度上有利于促進熔池內氣、液兩相間的交互作用。但同時也必須注意,熔體受氣流影響而獲得的運動速度越大,其對爐墻壁面可能造成的沖刷也更嚴重,甚至可能對金屬鉛的沉降造成影響,因此,在實際生產中選擇噴吹方式時,應當綜合考慮多方面因素的影響。

表5 不同高度以下空間區域內熔體的平均速度Table 5 Melt average velocity under different heights

2.2 不同噴吹速度對熔體運動的影響

2.2.1 工況參數設置

為研究不同一次風速對熔池內氣液兩相流動的影響,除生產中常用的225 m/s 的一次風風速外,還針對一次風速分別為175 m/s和275 m/s時的工況進行模擬計算。由于錯吹為目前生產中常用的噴吹方式,因此計算中采用上述錯吹工況的物理模型及物性參數條件。在錯吹方式下,不同風速工況參數設置如表6所示。

表6 錯吹方式下不同風速工況參數設置Table 6 Simulation parameters of different injection speeds

2.2.2 數值模擬結果分析

圖8所示為3個工況條件下高風口和低風口中軸線上的氣流貫穿深度隨時間的變化曲線。從圖8可以看出,在3個工況下,兩風口位置處的氣流貫穿深度隨時間變化的規律基本一致:在噴吹開始后的0~0.5 s 內,風口位置處的氣流貫穿深度波動均相對較大,但高風口位置處氣流貫穿深度的波動幅度略小于低風口位置處的結果,這是由于高風口在熔體中的浸沒深度較小,風口上部熔體所產生的靜壓也相對較小,氣流進入熔池后所產生的氣團溢出熔體表面的時間較短,且體積相對較小。

圖8 不同噴吹速度下氣流貫穿深度隨時間的變化曲線Fig.8 Change curves of air penetration depth with time at different injection speeds

從圖8還可以發現,無論是高風口處還是低風口處,氣流貫穿深度在1~2 s時均基本穩定在0.4 m左右。3個工況下1~2 s內氣流貫穿深度的最大值、最小值及平均值如表7所示。可見:無論在高風口還是在低風口處,氣流貫穿深度的最大值、最小值以及平均值均隨一次風速的提升而增大,且低風口處的氣流貫穿深度均大于高風口處的結果。這說明增大一次風速能夠在一定程度上提高氣流的貫穿深度,進而增大富氧空氣與熔渣之間的接觸面積,有利于反應的發生。

表7 不同工況下1~2 s內氣流貫穿深度的最大值、最小值及平均值Table 7 The maximum,minimum and average values of air penetration depth during 1-2 s氣流貫穿深度/m

圖9所示為3個工況下不同爐體高度位置截面上熔體相的平均體積分數的變化曲線。從圖9可以看出,在不同風速的條件下,熔體相體積分數的變化趨勢相似;但在靜止液面以下0.25 m范圍內,3個工況的熔體相體積分數均快速下降,這說明該區域內的大量熔體被氣流帶動翻卷至靜止液面以上的位置,即受一次風擾動較強烈的區域。另外,在該強烈擾動區域內,隨著風速增加,熔體相體積分數下降更快,這說明風速由175 m/s 提升至275 m/s 時,能夠為該區域內的熔體提供更強的擾動。

圖9 熔體相體積分數沿爐體高度方向上的變化曲線Fig.9 Melt volume fraction changing along height of furnace

圖10 所示為同風速下熔體區域的速度矢量圖(開放右側低風口)。從圖10可以發現,在一次風口上部,即使該截面位置為右側低風口噴吹,在左側仍會出現熔體向熔池中部翻卷的狀態,說明左側位置熔體會受到相鄰位置處一次風的帶動,出現熔體向熔池中部翻卷的現象。但是,隨著一次風速不斷增大,翻卷至爐體中部的熔體逐漸貼近左側爐墻,圖中紅色箭頭所示漩渦也越貼近左側爐墻,且當一次風速達到275 m/s 時一次風口下部的渦消失,即風速越高,右側(風口側)熔體受氣流帶動向左側翻卷的程度越劇烈,這說明風速越高,對熔池內熔體的影響區域越大。

圖10 熔體區域速度矢量圖Fig.10 Vector in melt zone

根據圖5(a)中所示不同高度位置,計算3 個工況中不同高度以下空間區域的熔體相平均速度,結果如圖11 所示。從圖11 可以看出,渣液面高度以下熔體的運動速度最大,低風口高度以下熔體運動速度其次,說明一次風進入熔池內部后對風口高度至渣液面高度空間區域的熔體運動影響最大。在相同高度位置下,熔體的運動速度隨著風速的增大而增大,說明在175 m/s至275 m/s的風速范圍內,增大風速能夠提高熔池內熔體運動的劇烈程度。另外,當風速由175 m/s 增大至225 m/s時,渣液面以下區域的熔體相平均速度變化劇烈,而鉛液面以下區域的熔體相的平均速度變化很小,說明增大風速對風口上部區域熔體的影響遠比對底部金屬相的影響更大,因此,在合理范圍內增大風速,能夠提高對一次風口上部熔體的擾動,從而提高氣流和熔體的接觸面積,有利于反應的發生和進行;同時,由于鉛液面以下金屬相的平均速度基本穩定不變,因此,風速的提高對于熔池底部金屬相的靜置不會產生嚴重的影響。

圖11 不同工況下不同高度位置以下空間區域內熔體的平均速度Fig.11 Melt average velocity under different heights in different cases

3 結論

1)在對吹工況下,低風口截面處熔體向熔池中部翻卷并基本呈對稱的狀態;在錯吹工況下,因受相鄰高風口帶動的影響在低風口截面處仍會出現與對吹時相似的結果。

2)在相同的風量條件下,采用4個低位風口對吹時氣流對熔體的影響更大,能夠使熔池中熔體運動更加劇烈,但熔體的運動速度越大,對爐墻壁面造成的沖刷蝕損可能更加嚴重,甚至可能影響金屬鉛的沉降過程,在實際生產中選擇噴吹方式時應當綜合考慮。

3)不同風速工況下,無論在高風口還是在低風口處,氣流貫穿深度隨風速的提升而增大,且低風口處的氣流貫穿深度均大于高風口處的結果。

4)在風速為175~275 m/s范圍內,一次風對熔體擾動強烈的區域為靜止渣液面至液面以下0.25 m范圍,且增大風速有利于提高熔池上部熔體運動的劇烈程度,進而提升熔煉效率,且不會對金屬相的靜置產生顯著影響。

猜你喜歡
風速區域
永久基本農田集中區域“禁廢”
今日農業(2021年9期)2021-11-26 07:41:24
分割區域
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
基于GARCH的短時風速預測方法
關于四色猜想
分區域
考慮風切和塔影效應的風力機風速模型
電測與儀表(2015年8期)2015-04-09 11:50:06
基于嚴重區域的多PCC點暫降頻次估計
電測與儀表(2015年5期)2015-04-09 11:30:52
GE在中國發布2.3-116低風速智能風機
主站蜘蛛池模板: 3344在线观看无码| 亚洲欧美国产视频| 欧美精品一区二区三区中文字幕| 国产AV毛片| 真实国产精品vr专区| 99re在线视频观看| 国产农村精品一级毛片视频| 香蕉eeww99国产在线观看| 91系列在线观看| 91精品国产91欠久久久久| 国产99欧美精品久久精品久久| 亚洲人成影院在线观看| 男女精品视频| 国产毛片高清一级国语| 精品成人一区二区三区电影| 色综合综合网| 成人日韩精品| 国产va欧美va在线观看| 色婷婷电影网| 99国产精品免费观看视频| 国产精品爆乳99久久| 国产极品嫩模在线观看91| 毛片基地美国正在播放亚洲| 四虎亚洲精品| 2021最新国产精品网站| 黄色在线网| 国产精品香蕉在线观看不卡| 婷婷在线网站| 欧美成人精品在线| 国产精品无码AV片在线观看播放| 国产精品美乳| 国产新AV天堂| 啪啪免费视频一区二区| 欧美日韩国产综合视频在线观看| 国产福利在线免费观看| 中文字幕乱码二三区免费| 欧美中文字幕一区| 中文字幕乱码二三区免费| 午夜一区二区三区| 国产欧美日韩专区发布| 91精品伊人久久大香线蕉| 2020最新国产精品视频| 九九精品在线观看| 2024av在线无码中文最新| 国产黑丝一区| 欧美日韩精品一区二区视频| 亚洲成人网在线观看| 亚洲美女视频一区| 国产精品不卡永久免费| 在线视频亚洲色图| 五月婷婷综合网| 一本久道久久综合多人| 影音先锋亚洲无码| 成人免费午夜视频| 午夜国产在线观看| 99久久精品国产麻豆婷婷| 伊人激情综合| 亚洲中文字幕97久久精品少妇| 中文字幕调教一区二区视频| 另类重口100页在线播放| 五月婷婷综合色| 国产美女在线观看| 亚洲一区无码在线| A级毛片高清免费视频就| 国产情精品嫩草影院88av| 女人18一级毛片免费观看| 色欲色欲久久综合网| 一区二区三区在线不卡免费| 亚洲精品另类| 在线免费观看a视频| 无码免费视频| 欧美成在线视频| 国产免费好大好硬视频| 久久婷婷国产综合尤物精品| 成人国产小视频| 天堂va亚洲va欧美va国产| 欧美成一级| 亚洲成网777777国产精品| 欧美一级色视频| 毛片免费试看| 国产日本欧美在线观看| 亚洲 日韩 激情 无码 中出|