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

固旋閥塔板的流體力學性能研究及其旋轉流場CFD模擬

2017-01-21 02:19:58王海鵬朱菊香姚克儉
石油化工 2016年9期

王海鵬,朱菊香,齊 亮,姚克儉

(浙江工業大學 化學工程學院 綠色化學合成技術國家重點實驗室培育基地,浙江 杭州 310032)

固旋閥塔板的流體力學性能研究及其旋轉流場CFD模擬

王海鵬,朱菊香,齊 亮,姚克儉

(浙江工業大學 化學工程學院 綠色化學合成技術國家重點實驗室培育基地,浙江 杭州 310032)

提出了一種新型固定旋轉閥塔板,在內徑600 mm的有機玻璃塔內,以空氣-水為物系,對固旋閥塔板的流體力學性能進行研究。測定了塔板壓降、漏液率和霧沫夾帶量等流體力學性能參數,并與旋轉浮閥塔板和F1浮閥塔板進行對比。實驗結果表明,當氣體負荷較大時,固旋閥塔板的干板壓降大于F1型浮閥塔板,具有更高的傳質效率;固旋閥塔板的濕板壓降小于旋轉浮閥塔板和F1浮閥塔板;固旋閥塔板的霧沫夾帶比F1型浮閥塔板小30%~40%,比旋轉浮閥塔板小10%~20%,具有更高的氣相負荷操作上限。通過Fluent6.3軟件的模擬計算分析,固旋閥塔板的流場更穩定,液層分布更均勻,操作性能更優。

固旋閥塔板;旋轉浮閥塔板;旋轉流場;壓降;漏液;霧沫夾帶

按照鼓泡元件的不同,板式塔主要可以分為泡罩塔、篩板塔、浮閥塔和固定閥塔等[1],其中,浮閥型塔板[2]由于氣體流通面積可隨氣體負荷的變化自動調節,具有良好的操作彈性,又由于改變了以往的塔板結構,減小了浮閥塔板上的液相返混,使氣流以水平方向吹向浮閥板面的液層[3-5],是一類高效的氣液傳質設備,在實際應用中占據著重要地位。

浙江工業大學對典型的F1浮閥塔板和F1浮閥全開后焊接在塔板上的傳質性能進行了研究,發現F1浮閥塔板的微旋轉運動可提高塔板傳質效率?;谏鲜鰧嶒灲Y果,浙江工業大學開發了一種新型傳質鼓泡元件——旋轉浮閥。張緒滿等[6]研究發現,旋轉浮閥塔板的霧沫夾帶率比F1型浮閥塔板低50%以上,壓降和F1型浮閥塔板壓降基本相同,漏液分率比F1型閥塔板約低21.19%;旋轉浮閥塔板的操作彈性遠大于F1型浮閥塔板。袁云峰等[7]研究結果表明,旋轉浮閥塔板的氣含率分布比F1型浮閥塔板均勻,且旋轉浮閥塔板的氣含率分布受液體噴淋密度和氣體動能因子的影響較小。但旋轉浮閥塔板仍存在著一些問題,一方面有浮閥旋轉易磨損,浮閥更易脫落等缺點[8-10],另一方面,浮閥在液層局部旋轉使得整個液層鼓泡不均勻,泡沫層高度分布不一,導致流場不穩定,最佳操作性能受到限制。

從優化塔板結構﹑強化旋轉流場﹑提高塔板效率等方面考慮,針對旋轉浮閥的不足之處,圍繞旋轉閥的概念,并結合固定閥的良好應用[11],本工作提出了一種新型的固定旋轉閥塔板,簡稱固旋閥塔板,對固旋閥塔板的流體力學性能進行考察研究,將性能測試結果與旋轉浮閥塔板和F1型浮閥塔板進行比較,探討氣液兩相流的新特性,并利用Flunet6.3軟件對固旋閥塔板的旋轉流場進行了計算流體動力學(CFD)單氣相模擬,探究板上氣相流場分布,從微觀上討論及分析固旋閥塔板的性能特征。

1 實驗部分

1.1 塔板結構與特點

圖1為固旋閥的結構。由圖1可知,該固旋閥是一種組合閥,由主閥部分和副閥部分組合而成。其中,主閥部分是圓形固定閥[12],其在塔板上的投影面為φ40 mm的圓形區域,副閥部分是1個閥蓋,其在塔板上的投影面為φ50 mm的圓形區域,閥蓋周邊設置向下彎曲的折邊,周圍均勻排布3個條形翅片。

固旋閥和旋轉浮閥一樣能夠在塔板液層內形成旋轉流場。當氣體從塔板下方穿過閥孔沖擊固旋閥時,由于閥蓋存在阻擋作用,使氣流從垂直運動變為水平運動。彎曲的折邊使氣體從閥的側孔斜向下吹到塔板板面上,可有效避免與鄰近閥孔吹出來的氣流發生直接對沖,有效降低霧沫夾帶和塔板漏液,提高塔板的處理能力和操作彈性。

翅片具有較強的導向作用,可使氣體從閥的側孔以螺旋旋轉的方式流出,并機械性地分割氣液兩相,使氣液更為細化,大幅增加了氣液的接觸面積。同時,導向旋轉運動可以加強塔板上的氣液湍流程度,加快閥體附近氣液表面更新,使得氣液接觸更為均勻并減少死區,提高傳質效率。

圖1 固旋閥的結構Fig.1 Structure of a fxed rotary valve.

1.2 實驗裝置

在透明的有機玻璃塔中采用空氣-水體系進行冷模實驗,測試了固定旋轉閥塔板的流體力學性能,實驗裝置見圖2。塔內安裝3塊相同的塔板,塔板的結構參數見表1。其中,中間層塔板為測試塔板,下層塔板為氣體分布板,上層塔板為液體分布板,用弓形降液管作為降液通道。此外在上層塔板之上安裝霧沫夾帶收集板,用于收集霧沫夾帶量,同時塔的頂部安裝一層絲網填料除霧器,以保證準確測量霧沫夾帶量;在下層塔板之下安裝一層漏液收集板兼做氣體初步分布板。

表1 實驗塔及塔板的結構參數Table 1 Structure parameters of the experimental column and trays

實驗過程中,空氣由離心式鼓風機提供,經對夾式氣體孔板流量計后從塔底進入塔內,氣量由U形壓差計讀數換算得到;水由離心泵從塔釜輸送至塔頂,最底部塔板流出的液體返回塔釜循環使用,液量由轉子流量計測定。塔板壓降采用U形壓差計測量,塔板的漏液量和霧沫夾帶量采用液體收集器進行收集并用電子秤計量。

圖2 實驗裝置Fig.2 Schematic diagram of the experimental column.1 Centrifugal blower;2 Orifce-plate fowmeter;3,4,12 Diferential manometers;5 Collecting tray of weeping;6,8 Trays;7 Experimental tray;9 Collecting tray of entrainment;10 Wire-mesh demister;11 Pump

2 結果與討論

2.1 干板壓降

干板壓降是氣體通過塔板時遇到阻力而產生的能量損失。固旋閥塔板,旋轉浮閥塔板和F1浮閥塔板的干板壓降(Δpd)與閥孔動能因子(F0)的關系見圖3。

圖3 固旋閥塔板、旋轉浮閥塔板和F1浮閥塔板的干板壓降比較Fig.3 Comparison of dry pressure drops of the fxed rotary valve tray,rotary foat valve tray and F1 valve tray.Δpd:dry pressure drop;F0:the kinetic energy factor of valve hole.■ Fixed rotary valve tray;● Rotary foat valve tray;▲ F1 valve tray

由圖3可知,3種塔板的干板壓降均隨著F0的增大而逐漸增加。其中,旋轉浮閥塔板和F1型浮閥塔板均存在開啟前、開啟中和全開后3個階段。當F0>17 (m·s-1)(kg·m-3)0.5時,固旋閥塔板的干板壓降大于F1浮閥塔板,但比旋轉浮閥塔板小,這是由于固旋閥沒有像旋轉浮閥那樣的浮動部件,減少了用于抬起浮閥的能量損失;而相對于F1浮閥塔板,由于固旋閥翅片的存在而形成了旋轉流場,氣體通過閥孔后流向發生了變化,損失更多動能,表現為其干板壓降稍大于F1浮閥塔板。

采用Prince關聯式(1)[13]56-65:

式中,Δpd為干板壓降,Pa;ρG,ρL分別為氣相和液相的密度,kg/m3;u0為閥孔氣速,m/s;C0為孔流系數。

對固旋閥塔板以及全開后的旋轉浮閥塔板和F1型浮閥塔板進行干板壓降關聯,得到各塔板的孔流系數值,按照關聯式(1)計算出的各塔板干板壓降計算值和實驗值進行對比,計算值和實驗值吻合較好。表2為各塔板的孔流系數和決定系數。由表2可知,孔流系數越小,其閥孔的流動阻力越大,產生的干板壓降就越大,與實驗結果相一致。

表2 各塔板的孔流系數和決定系數Table 2 Orifce coefcients andR2of the diferent trays

2.2 濕板壓降

濕板壓降是指氣相通過塔板和板上清液層時產生的阻力損失,是流體力學性能的重要參數之一。本工作主要研究了氣速和噴淋密度(L)變化對實驗塔板濕板壓降的影響,圖4為固旋閥塔板的濕板壓降(Δpw)和F0的關系。由圖4可知,在F0相同的條件下,Δpw隨L的增加而增大,這是因為塔板上的清液層高度不斷增加,氣體通過液層的阻力也相應逐漸增大;當F0>22(m·s-1)(kg·m-3)0.5時,各L下的Δpw數值差距越來越小。這主要是由于氣速很大時,塔板上液體均被大量吹起,各L的塔板清液層高度接近,此時的Δpw主要是由Δpd和氣體通過厚度差距不大的液層產生的壓降組成。結合實驗現象可知,當F0較大時,不同的L下,板上的清液層高度較為相近。

圖4 固旋閥塔板的Δpw與F0的關系Fig.4 Relationship between ΔpwandF0of the fxed rotary valve tray.L(spray density)/(m3·m-2·h-1):■ 20;● 30;▲ 40;▼ 60 Δpw:wet pressure drop.

計算濕板壓降的關聯式主要有加和式、準數關聯式和氣速關聯式。通常情況下,準數關聯式更能清晰表明其關鍵影響因素。對于空氣-水體系,計算濕板壓降的關聯式主要采用式(2)[13]357-360。

式中,Δpw是濕板壓降,Pa;α為系數;β1,β2,β3為指數;L為噴淋密度,m3/(m2·h);F0為閥孔動能因子,(m·s-1)(kg·m-3)0.5;hw為溢流堰高度,m;由于實驗過程中hw是一個固定值,所以式(2)可簡化為式(3)。

采用關聯式(3),用最小二乘法對Δpw數據進行擬合,得到回歸方程,見式(4)~(5)。

式中,F0的適用范圍為:4.42~11.76(m·s-1)·(kg·m-3)0.5,決定系數R2=0.989。

式中,F0的適用范圍為:11.76~25.62(m·s-1)·(kg·m-3)0.5,決定系數R2=0.991。

圖5為固旋閥塔板、旋轉浮閥塔板和F1浮閥塔板的Δpw比較。由圖5可知,固旋閥塔板的Δpw小于旋轉浮閥塔板和F1浮閥塔板。當F0<11(m·s-1)·(kg·m-3)0.5時,固旋閥塔板的Δpw比旋轉浮閥塔板和F1浮閥塔板小,此時板上清液層高度相近,板壓降的不同主要是由于Δpd的不同而造成的;當F0>11(m·s-1)(kg·m-3)0.5時,固旋閥塔板Δpw比F1浮閥塔板小,一方面是由于固旋閥沒有活動部件,少了抬起浮閥的能量損失,另一方面,由于翅片具有較強的導向作用,使氣體從閥的側孔以螺旋旋轉的方式流出,并機械性地分割氣液兩相,使氣液兩相更為細化,固旋閥塔板的Δpw小于F1型浮閥塔板,同時,由于固旋閥塔板具有更為穩定的旋轉流場,導向旋轉運動加強塔板上的氣液湍流程度,加快閥體附近氣液表面更新,傳質鼓泡分布更為均勻,降低了清液層高度,因此,固旋閥塔板的Δpw也小于旋轉浮閥塔板。

圖5 固旋閥塔板、旋轉浮閥塔板和F1浮閥塔板的Δpw比較Fig.5 Comparison of Δpwof the fxed rotary valve tray,rotary foat valve tray and F1 valve tray.■ Fixed rotary valve tray;● Rotary foat valve tray;▲ F1 valve trayL/(m3·m-2·h-1)=40.

2.3 霧沫夾帶

霧沫夾帶是指氣體負荷較高時下一層塔板上的液體以液滴的形式被氣體吹到上一層塔板,形成返混的過程。圖6為固旋閥塔板的霧沫夾帶(ev)和F0的關系。由圖6可知,在相同L下,固旋閥塔板的霧沫夾帶隨著F0的增大而逐漸增大。當F0相同時,隨著L的增加固旋閥塔板的ev逐漸增大;當F0較小時,固旋閥塔板的ev趨向于零,這也符合固閥類塔板ev的規律;當F0較大時,固旋閥塔板的ev急劇增加。由圖6還可知,當L較大時,固旋閥塔板的ev隨F0的增加而快速增加。

圖6 固旋閥塔板的ev和F0的關系Fig.6 Relationship between entrainment(ev) andF0of the fxed rotary valve tray.L/(m3·m-2·h-1):■ 20;● 30;▲ 40;▼ 60

圖7為固旋閥塔板、旋轉浮閥塔板和F1型浮閥塔板ev的比較。由圖7可知,在L相同時,固旋閥塔板的ev低于旋轉浮閥塔板和F1浮閥塔板。其中,當F0達到27(m·s-1)(kg·m-3)0.5時,旋轉浮閥塔板的ev達到上限,F1型浮閥塔板早已霧沫夾帶液泛,而固旋閥塔板的ev僅為0.1左右。這主要是由于彎曲的折邊使氣體從閥的側孔斜向下吹回到塔板板面上,可避免鄰近閥孔吹出來的氣流發生直接對沖,又存在翅片的導流作用,更好的消除液面梯度,使鼓泡更為均勻,液滴尺寸更為細小,并且氣體圍繞翅片發生旋轉作用,增長了氣體的流道,減小了氣體的動能,從而有效降低地霧沫夾帶,使固旋閥塔板的ev遠小于F1浮閥塔板,小30%~40%。另一方面,由于固旋閥良好的結構設計,強化了旋轉流場,使液層分布更為均勻,達到均一化。結合實驗現象可知,當F0較大時,固旋閥塔板的泡沫層高度低于旋轉浮閥塔板,所以固旋閥塔板的ev小于旋轉浮閥塔板,小10%~20%。

圖7 固旋閥塔板、旋轉浮閥和F1浮閥塔板的霧沫夾帶的比較Fig.7 Comparison of the entrainments of the fxed rotary valve tray,rotary foat valve tray and F1 valve tray.■ Fixed rotary valve tray;● Rotary foat valve tray;▲ F1 valve trayL/(m3·m-2·h-1)=40.

2.4 塔板漏液

塔板泄漏是指當閥孔氣速小于某一定值時,氣體不能支撐住塔板上的液體,使液體從閥孔大量漏至下一層塔板,造成液相返混,塔板效率下降的現象。圖8為固旋閥塔板的漏液分率和F0的關系。由圖8可知,在L相同時,固旋閥塔板的漏液分率隨著F0的增大而減?。辉贔0相同時,隨著L的增加,固旋閥塔板的漏液分率不斷減小,這一規律與其他固閥類塔板的漏液規律相同。當F0<6.5 (m·s-1) ·(kg·m-3)0.5時,固旋閥塔板的漏液分率隨著F0的增大而急劇減小,各L的漏液分率線都很清晰的分開;但當F0>6.5(m·s-1)(kg·m-3)0.5時,不同L下的漏液分率線逐漸趨向一致或相互交叉,這是因為當F0增大時,固旋閥塔板的漏液量逐漸減小,所以不同L下的漏液分率非常接近。

圖8 固旋閥漏液分率與F0的關系Fig.8 Relationship between the weeping fraction andF0ofthe fxed rotary valve tray.L/(m3·m-2·h-1):■ 20;● 30;▲ 40;▼ 60

圖9 為固旋閥塔板、旋轉浮閥塔板和F1浮閥塔板的漏液分率的比較。由圖9可知,在L相同時,固旋閥塔板的漏液分率大于旋轉浮閥塔板和F1浮閥塔板,這是因為固旋閥塔板的實際開孔率不會像浮閥塔板(沒有自動調節的功能)那樣隨著氣流的變化而發生變化。因此,當R2較小時,氣體不能完全托住液體,此時漏液量較大。

圖9 固旋閥塔板、旋轉浮閥塔板和F1浮閥塔板的漏液分率的比較Fig.9 Comparison of weeping fractions of the fxed rotary valve tray,rotary foat valve tray and F1 valve tray.■ Fixed rotary valve tray;● Rotary foat valve tray;▲ F1 valve trayL/(m3·m-2·h-1)=40.

3 旋轉流場CFD模擬

為了進一步探究固旋閥塔板的旋轉流場對其流體力學性能的影響,通過Fluent6.3軟件對固旋閥塔板上氣相三維流場進行了數值模擬,重點研究了塔板上固旋閥周圍以及相鄰閥件之間的速度分布,更系統的了解固旋閥塔板上的氣相流動特征。采用Gambit按照實際尺寸和形狀建立固旋閥物理模型,靠近閥件及塔板處采用非結構四面體網格進行劃分,其他區域采用六面體網格。Fluent6.3模擬計算過程中采用壓力基隱式求解器;采用湍動黏度計算公式中引入了旋轉與曲率的Realizableκ-ε湍流模型[14];采用有限體積法離散控制方程;采用二階精度的Upwind離散格式處理對流項;采用SIMPLEC算法處理壓力—速度耦合;亞松馳因子使用Flunet6.3軟件中的默認值,收斂條件為10-5。入口邊界條件采用速度進口,假設入口處速度分布均勻,且與入口平面垂直,入口處的湍流參數用湍流強度0.16(Re)-1/8和水力學直徑描述;出口邊界條件采用壓力出口,假設出口方向的各個流動變量的擴散通量為0;壁面條件采用無滑移邊壁,在近壁面采用標準壁面函數。

圖10為固旋閥塔板干板壓降實驗值和模擬值的對比結果,由圖10可知,模擬值出現了一定的絕對偏差,但其相對誤差小于3.5%,可認為和實驗值吻合,模擬結果具有真實可靠性。

圖11為固旋閥塔板和F1浮閥塔板X-Y剖面(板上8 mm)上的氣相流場分布圖。由圖11可知,固旋閥上的翅片具有較強的導向作用,使氣體從閥的側孔以螺旋旋轉的方式流出,一方面,在閥件之間形成一個旋轉流場,加強了塔板上的氣液湍流程度,使得氣液接觸更為充分并減少死區,提高了傳質效率;另一方面,翅片改變了氣體的流動方向,有效地避免鄰近閥孔吹出來的氣流發生直接對沖,從而大幅降低了固旋閥塔板的霧沫夾帶。同時,由于旋轉浮閥處于開啟階段時可以自由地上下浮動進行微轉動,導致相鄰的閥件上翅片排布具有隨機性,故形成的流場分布不穩定,此處不做模擬計算,而固旋閥固定的結構特征使得臨近閥件的翅片可進行合理排布,強化了旋轉流場的作用,使其具有比旋轉浮閥更優的操作性能??梢?,微觀的模擬計算結果能夠較好地解釋宏觀的實驗測試結果。

圖10 固旋閥塔板的Δpd實驗值與模擬值的對比Fig.10 Comparison between measured values and simulated values for Δpdof the fxed rotary valve tray.■ Measured value;● Simulated value

圖11 塔板X-Y剖面(Z=8 mm)氣相流場分布圖Fig.11 Gas fow feld distributions of theX-Ysections(Z=8 mm) of the valve trays.(a) Fixed rotary valve tray;(b) F1 valve tray

4 結論

1)在傳統固定閥的基礎上,結合旋轉閥的理念,提出了一種結構新穎的固旋閥塔板。

2)當氣體負荷較大時,固旋閥塔板的干板壓降大于F1型浮閥塔板,具有更高的傳質效率;固旋閥塔板的濕板壓降小于旋轉浮閥塔板和F1浮閥塔板;固旋閥塔板的霧沫夾帶比F1型浮閥塔板小30%~40%,比旋轉浮閥塔板小10%~20%,具有更高的氣相負荷操作上限。

3)通過Fluent6.3軟件的模擬計算分析,Realizableκ-ε湍流模型可用于固旋閥塔板上旋轉流場的模擬與分析,其結果從微觀上解釋了固旋閥塔板的流體力學性能特征,固旋閥塔板的流場更穩定,液層分布更均勻,操作性能更優。

4)固旋閥塔板是一種大通量、高效、高彈性

的新型塔板,具有廣泛的應用前景。

符 號 說 明

[1]王樹楹,高長寶,蘭仁水,等. 板式塔研究進展[J]. 化學工程,2003,31(3):20 - 26.

[2]梁治國. BTV4浮閥塔板的開發研制[D]. 東營:中國石油大學(華東),2006.

[3]潘忠濱,姚克儉,李育敏,等. 新型浮閥塔板研究進展[J].化工進展,2005,24(9):956 - 963.

[4]王少鋒,項曙光. 浮閥塔板最新應用研究進展[J]. 化工進展,2014,33(7):1677 - 1683.

[5]李玉安,路秀林,趙培,等. 導向浮閥塔板和F1型浮閥塔板的比較[J]. 石油化工,1996,25(8):563 - 568.

[6]張緒滿,姚克儉,何建烽,等. 旋轉浮閥塔板的流體力學性能研究[J]. 石油化工,2012,41(8):916 - 920.

[7]袁云峰,齊亮,姚克儉,等. 旋轉浮閥塔板的氣含率分布研究[J]. 石油化工,2013,42(9):990 - 995.

[8]張杰旭,李玉安,趙培,等. B型導向浮閥塔板實驗研究[J].化學工程,2005,33(6):1 - 3.

[9]左美蘭. 新型浮閥塔板水力學性能的研究[D]. 青島:青島科技大學,2009.

[10]任敏山,姚克儉,李育敏,等. 齒邊浮閥塔板流體力學性能的研究[J]. 石油化工,2005,34(4):356 - 359.

[11]Nutter D E. The MVGTMtray at FRI[J]. Trans IChemE,A,1999,77(6):493 - 497.

[12]付小蘇,張頌紅,姚克儉,等. 圓形固定閥塔板的性能研究及其工業應用[J]. 石油化工,2011,40(7):775 - 779.

[13]蘭州石油機械研究所. 現代塔器技術[M]. 北京:中國石化出版社,2005:56 - 65,357 - 360.

[14]Shih T H,Liou W W,Shabbir A,et al. A newκ-εeddy viscosity model for high reynolds number turbulence flows[J]. Comput Fluids,1995,24(3):227 - 238.

(編輯 楊天予)

Hydrodynamic performances of fixed rotary valve tray and CFD simulation of rotary flow field

Wang Haipeng,Zhu Juxiang,Qi Liang,Yao Kejian
(College of Chemical Engineering,Zhejiang University of Technology,State Key Laboratory Breeding Base of Green Chemistry-Synthesis Technology,Hangzhou Zhejiang 310032,China)

A kind of novel fixed rotary valve trays was presented through several improvements of traditional fxed valve trays. The hydrodynamic performances of the fxed rotary valve trays were experimentally investigated with an air-water system in an organic glass column with an inner diameter of 600 mm. The pressure drop,weeping fraction and entrainment of the fxed rotary valve trays were measured and compared with those of rotary foat valve trays and F1 valve trays. The results showed that the dry pressure drop of the fxed rotary valve trays was higher than that of the F1 valve trays while the kinetic energy factor got larger. The wet pressure drop of the fxed rotary valve trays was lower than those of the rotary foat valve trays and F1valve trays. The entrainment of the fxed rotary valve trays was less than that of the F1 valve trays by 30%-40% and less than that of the rotary foat valve trays by 10%-20%. According to the results of the CFD simulation,the performances of the fxed rotary valve trays were much better than those of the rotary foat valve trays and F1 valve trays.

fxed rotary valve trays;rotary foat valve trays;rotary fow feld;pressure drop;weeping;entrainment

1000 - 8144(2016)09 - 1100 - 07

TQ 053.5

A

10.3969/j.issn.1000-8144.2016.09.013

2016 - 02 - 16;[修改稿日期]2016 - 06 - 03。

王海鵬(1990—),男,浙江省臺州市人,碩士生。聯系人:姚克儉,電話 0571 - 88320952,電郵 yaokj@zjut.edu.cn。

主站蜘蛛池模板: 粗大猛烈进出高潮视频无码| 精品国产99久久| 亚洲区欧美区| 粉嫩国产白浆在线观看| 国产黄色视频综合| 国产精品妖精视频| 黄片一区二区三区| 99草精品视频| 亚洲精品国产自在现线最新| 99精品视频播放| 午夜福利亚洲精品| 99久久这里只精品麻豆| 亚洲成在线观看| 欧美日本视频在线观看| 久久婷婷六月| 99精品福利视频| 亚洲区一区| 亚洲精品成人7777在线观看| 国产精品九九视频| 国产精品永久免费嫩草研究院| 欧美一区二区自偷自拍视频| 久久久久久久97| 亚洲中文字幕无码爆乳| 免费A级毛片无码无遮挡| 在线高清亚洲精品二区| 女人爽到高潮免费视频大全| 精品国产香蕉在线播出| 国产精品极品美女自在线| 欧美精品在线免费| 全色黄大色大片免费久久老太| 国产手机在线小视频免费观看 | 免费看a毛片| 婷婷综合缴情亚洲五月伊| 国产区免费| 免费国产高清视频| 五月婷婷亚洲综合| 色妞永久免费视频| 欧美在线黄| vvvv98国产成人综合青青| 欧美一级视频免费| 蜜桃视频一区二区| 欧美日韩第二页| 四虎精品黑人视频| 欧美一道本| 国产亚洲精| 狠狠干欧美| 午夜天堂视频| 欧洲高清无码在线| 国产在线日本| 全部免费毛片免费播放| 久久精品亚洲热综合一区二区| 人妻无码中文字幕第一区| 亚洲水蜜桃久久综合网站 | 九九久久精品免费观看| 国产成人91精品| a毛片免费看| 国产午夜无码专区喷水| 五月婷婷导航| 久久99久久无码毛片一区二区 | a级毛片免费播放| 亚洲色图欧美一区| 538精品在线观看| 亚洲一区二区日韩欧美gif| 91无码国产视频| 在线欧美a| 91亚洲精品第一| 青青草原国产免费av观看| 日韩精品亚洲精品第一页| 高清不卡一区二区三区香蕉| 中文字幕中文字字幕码一二区| 91精品网站| 国产日韩欧美在线视频免费观看| 亚洲自偷自拍另类小说| 72种姿势欧美久久久大黄蕉| 日韩免费毛片视频| 亚洲精品无码AⅤ片青青在线观看| 都市激情亚洲综合久久| 无码中文字幕加勒比高清| 一本综合久久| 日韩人妻少妇一区二区| 精品一区二区三区波多野结衣| 99久久精品免费观看国产|