• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    規(guī)整填料表面液體分布的計算流體力學

    2015-11-26 02:43:06譚麗媛袁希鋼KALBASSIMohammadAli
    化工進展 2015年9期
    關鍵詞:規(guī)整液膜表面張力

    譚麗媛,袁希鋼,KALBASSI Mohammad Ali

    ?

    規(guī)整填料表面液體分布的計算流體力學

    譚麗媛1,袁希鋼1,KALBASSI Mohammad Ali2

    (1天津大學化學工程聯(lián)合國家重點實驗室,天津 300072;2Air Products Public Limited Company,Walton on Thames,UK)

    采用計算流體力學(computational fluid dynamics,CFD)中的VOF方法對規(guī)整填料表面液相分布進行了三維建模和仿真,實現(xiàn)規(guī)整填料內液相分布的可視化并得到了液膜厚度和有效相界面積比等相關定量信息。通過分析比較不同物系的數(shù)值模擬結果,發(fā)現(xiàn)液體的表面張力和黏度都對填料表面上液體分布有影響。表面張力越小,液相在填料片上分布越均勻,有效相界面積比越大,液膜厚度越小;黏度增加,有效相界面積比和液膜厚度也隨之增加。較之于黏度對液相分布的影響程度,表面張力的影響程度更大,為主要影響因素。本文還提出一個預測有效相界面積比的新公式,并將數(shù)值模擬結果與已有文獻進行對比,吻合性較好。

    計算流體力學;可視化;液膜厚度;有效相界面積比;表面張力;黏度

    規(guī)整填料塔因在汽提和吸收過程中為氣液兩相提供有效的接觸面而被廣泛應用于工業(yè)生產中,其分離效率的預測和提高一直以來都是研究人員最為關注的問題。Mohamed等[1]在對氣液相在規(guī)整填料塔內分布的研究中發(fā)現(xiàn),液體在規(guī)整填料塔內的不均勻性分布是影響塔分離效率的決定性因素,而氣相即使是在初始時就分布不均勻,也不會對塔的分離效率有影響。因此規(guī)整填料表面液相分布的研究在規(guī)整填料基礎研究中占有重要地位。液體在小尺度區(qū)域(基本結構單元)內的流動會在很大程度上影響液體在宏觀尺度上的分布。在一個結構單元內,液體的分布是由液體流過填料波紋波峰、填料節(jié)點和填料表面的孔洞所決定的。然而對結構單元內液體流動分布的定量特征和平均液膜厚度的研究卻鮮有研究報道[2]。

    早期對規(guī)整填料內液體宏觀流動分布的研究模型多為基于實驗研究提出的經驗或半經驗理論模型,如針對板波紋填料提出的SRP模型[3-6]、基于規(guī)整填料波紋幾何形狀Olujic等[7-8]提出的Delft模型、徐崇嗣等[9]通過對Dangizer[10]提出的靜態(tài)混合器中液體分布模型進行分析研究而建立的節(jié)點網格模型。針對規(guī)整填料內液體的微觀流動,如規(guī)整填料表面液膜流動,Nusselt[11]是最早對液膜流動進行研究的學者之一,他提出的降膜理論到現(xiàn)在仍可用于估算規(guī)整填料內的平均液膜厚度。Shetty和Cerro[12-13]最早提出了用于規(guī)整填料塔設計的一套完整的液膜流動模型。雖然液膜流動模型是以表面帶有波紋的規(guī)整填料為研究背景,但與實際填料上的液膜流動情況仍有很大差距。谷芳[14]、陳江波[15]采用LDV技術對規(guī)整填料表面液相流動分布進行了測試分析。

    最近幾年,隨著計算流體力學(computational fluid dynamics,CFD)和相關學科的發(fā)展及其能有效地簡化實驗設備和減少相關設備成本,CFD已成為研究復雜模型中多相流必不可少的工具[16-17]。另外,采用實驗的方法不能對填料塔內的流體力學參數(shù)提供詳細的信息,也很難利用探針測量規(guī)整填料塔內下降液膜的速度、溫度和壓力分布,這也使得數(shù)值模擬在這方面更加重要[18]。

    受目前計算能力的限制,對整個填料塔內每一現(xiàn)象(流體動力學、傳質、化學反應)進行模擬是非常困難的。因此,對規(guī)整填料塔的數(shù)值模擬研究主要從3個不同的尺寸:微尺度、中尺度和宏觀尺度[19-20]進行研究分析。微尺度的模擬通常是毫米或者厘米尺寸的,計算域網格大小是微米級的。Szulczewska等[21]基于Mellapak 250.Y型板波紋規(guī)整填料建立一個二維兩相流CFD模型,模擬了規(guī)整填料簡化模型上氣液兩相逆流流動過程。二維模擬不僅可以大大簡化數(shù)值模擬方法和顯著減少計算量,當液相以完全膜狀流進行流動時,二維模擬還可以對其流體動力學參數(shù)給出很好的近似。但隨著液體流率的減小,表面張力的作用會使液體流動由溪流變?yōu)榈瘟?,這些行為均為三維結構,無法使用二維模擬進行模擬觀察[22]。Iso等[23-24]對傾斜板上液膜流動進行了三維CFD模擬,很好地觀察了液體的溪流和滴流等行為。對于微尺度,液體的不均勻分布流體動力學是影響規(guī)整填料塔性能的關鍵因素,對其進行研究是很有必要的。中尺度的研究包含氣、液相在金屬規(guī)整填料內的流動,可以給出填料塔的干板壓降等流體力學參數(shù)。在研究過程中,可以用規(guī)整填料一系列的特征單元代表整個填料。Hosseini等[25]利用不同的湍流模型研究了氣液兩相在Mellapak Plus 752.Y規(guī)整填料上的流體力學特性。Haroun等[26]對具有不同接觸角的液體在Mellapak 250.Y規(guī)整填料特征單元上液相流動分布情況進行了三維數(shù)值模擬,并對其有效相界面積和持液量進行了預測。Shojaee等[27]利用VOF模型研究了Gempak2A規(guī)整填料特征單元上液相流動分布及氣、液相流率對有效相界面積的影響。Dai等[28-29]采用CFD中歐拉多相流模型對規(guī)整催化填料結構與其壓降和傳質系數(shù)的關系進行了三維數(shù)值模擬,并通過優(yōu)化規(guī)整催化填料的結構參數(shù)提高其傳質系數(shù)和降低壓降;之后又利用此模型對泡點反應器中苯與丙烯烷基化反應進行了模擬,確定了規(guī)整催化填料的最佳操作條件。宏觀尺度利用從中尺度研究得到的相關關聯(lián)式作為輸入,進而決定填料塔內的流型以及壁面、分布器等對其影響。Raynal等[19]運用從全規(guī)模反應中尺度下得到的壓降關聯(lián)式檢查3種不同幾何構型的速度模型:無氣相分布器、垂直管分布器、沖擊板垂直分布器;最后用多孔介質模型模擬全塔,并得到較好的結果。

    綜上所述,早期對規(guī)整填料塔內流動的研究主要集中在宏觀流動分布方面,目前則集中在填料內部及更小尺度的流動描述,而且采用CFD方法研究填料內部的流動成為一種趨勢。雖然已有文獻針對不同尺度的液相不均勻分布的流體動力學進行研究,但是這兩者之間的特點仍需要在微尺度方面進行全面的理解。本文致力于研究液相物性對液體在填料表面分布的影響,建立了基于VOF方法的三維CFD數(shù)值分析模型。在研究過程中,作者首先對不同表面張力和黏度的物系在規(guī)整填料簡化模型上的分布進行分析研究,在此基礎之上,對液體在Mellapak 500.Y規(guī)整填料的近似結構單元上的分布進行模擬研究。

    1 數(shù)值模擬模型

    1.1 物理模型

    基于規(guī)整填料的結構特點,傾斜平板(或波紋板)被認為是規(guī)整填料片的最簡化單元。研究液體表面張力和黏度對液體在傾斜板上的分布是后續(xù)研究液體在規(guī)整填料更復雜模型上分布的一個重要步驟。

    本文所采用的傾斜波紋板的計算模型如圖1所示,該物理模型和流動條件的設置是基于現(xiàn)有傾斜板上液體流動分布研究[22-24]來設定的。物理模型尺寸為60mm×50mm×7mm,板上小波紋間距= 2.8mm,波峰=0.3mm,波紋傾角=60°。由于液相雷諾數(shù)較低(<300),液體在波紋板上流動方式視為層流,利用Nusselt液膜理論[11]估算液相進口區(qū)域為50mm×0.4mm。模擬計算過程中液相進口邊界條件取為速度進口 Velocity-inlet,L=0.3m/s;左右兩邊側壁和波紋板為不銹鋼材質且邊界條件取無滑移壁面;出口邊界條件取為壓力出口Pressure- outlet,出口壓力為0。為確定波動液膜的自由表面位置,采用結構化網格劃分該模型,并對近壁面區(qū)域進行網格加密,同時選用帶有幾何重建格式的VOF模型進行相界面的追蹤和捕獲。

    精餾塔規(guī)整填料通常為波紋板或者絲網結構。在填料塔中,規(guī)整填料被排成較塔內截面略小的圓盤[圖2(a)]),一盤之內,左右相鄰兩片規(guī)整填料的波紋傾斜方向相反。金屬板波紋填料上波紋與長邊的斜角為45°或60°,板片上打出一些小孔以利于氣體穿過[圖2(b)]。為減少計算時間,對規(guī)整填料結構單元進行一定的簡化,規(guī)整填料片表面上沒有開孔,三維CFD數(shù)值模型如圖2(c)所示。模型上方即(1)設為液相速度進口(velocity-inlet)邊界條件,L=0.01m/s;模型下方即(2)設為壓力出口(pressure-outlet)邊界條件,出口壓力為0;前后兩片填料為不銹鋼材質即(3)設為無滑移壁面 (wall);為了保證液體流動的連續(xù)性,模型左右兩邊即(4)有一定的液相流入,可近似為周期性邊界條件。規(guī)整填料的結構尺寸=10.10mm,波高=6.52mm。由于規(guī)整填料物理模型結構較為復雜,采用非結構化網格對其進行劃分,近壁面處進行網格加密,選用帶有幾何重建格式的VOF模型進行氣液相界面的追蹤和捕獲。

    1.2 數(shù)學模型

    本文采用的數(shù)值模擬方法是基于Hirt等[30]提出的適用于追蹤兩非混相界面的VOF方法。在流體動力學模型中,利用數(shù)值求解質量守恒和動量守恒。質量守恒可導出連續(xù)性方程,動量守恒可導出Navier-Stokes方程,即連續(xù)性方程如式(1)。

    Navier-Stokes方程如式(2)。

    (2)

    式中、、、和為局部的速度、密度、壓力、黏度和重力加速度;為除壓力、黏性力和重力之外的作用于流體上的其他力,本文中取表面張力。在研究過程中,采用Brackbill等[31]提出的連續(xù)表面力(continuum surface force,CSF)模型,使得表面張力在動量方程中以源相的形式加入到VOF的計算中,如式(3)。

    式中,l為網格單元內的液相體積分數(shù),0

    (4)

    為表面張力系數(shù),為自由表面曲率,其定義式如式(5)。

    (6)

    式中,、w為壁面處的單位法向量和切向量;w為壁面處界面切線與壁面的夾角。液相與壁面之間的接觸角及液相的表面張力利用OCA15EC光學接觸角測量儀測量得到。

    無論是對氣相還是液相,每一相的體積分數(shù)都遵從該相的轉運方程,如式(7)。

    (8)

    1.3 數(shù)值模擬策略

    本文主要考察液相在規(guī)整表面及其簡化模型上的分布及液體物性對其分布的影響情況。為了簡化數(shù)值模擬計算的時間,在模擬過程中假定氣相在填料內沒有流動,液體的流動方式為層流。即氣體速度為零,對液體的流動沒有影響。已有的研究表 明[14,27],載點以下氣相對液相的影響作用很小,因此在模擬過程中假定氣相為停滯相是合理的。為模擬填料表面瞬態(tài)氣液兩相流體動力學特性,VOF模型采用顯示差分格式,采用界面重構方法(geo-reconstruct method)作為體積分數(shù)離散格式,二階迎風格式求解動量方程,壓力項選用PRESTO!算法,壓力-速度耦合方程的求算采用PISO格式。計算過程中,壓力項的松弛因子為0.3,動量項為0.45,若連續(xù)性方程的殘差降至1×10-4以下,各個方向速度的殘差降至1×10-5時,可認為方程已收斂。模型計算推進時間步長Δ取1×10-5s,確保能準確模擬研究工況速度范圍內氣液兩相流體動力學特性。壁面為無滑移壁面,給定液相與壁面之間的接觸角。模擬計算初始時氣相充滿整個模型,即氣相的體積分數(shù)g=1,l=0。經網格獨立性驗證后,傾斜波紋板的網格數(shù)為233840,規(guī)整填料片的網格數(shù)為638935。

    1.4 模擬物系

    本文主要考察液體物性如黏度、表面張力等對液體在波紋規(guī)整填料片上的分布及潤濕情況的影響。表1列出了模擬物系的物理性質。

    表1 模擬物系的物性

    2 模型驗證

    為了驗證該計算模型和方法的準確性,按照文獻[23-24]中的模擬條件對進行了相應的模擬計算,對比了在相同情況下液相在傾斜波紋板面上的流動分布情況。圖3為兩種模擬結果的對比,二者吻合很好,證明了該數(shù)值模型和計算方法的準確性。

    3 液相在傾斜波紋板上的分布

    3.1 液相物性對液相流型的影響

    圖4給出了以平行于波紋板的視角液體在傾斜波紋板上的分布情況,顏色的深淺表征液膜厚度的大小。以液體在傾斜波紋板上的潤濕面積比率來考察液相在傾斜波紋板上的潤濕情況,潤濕面積比率為液相在板面上的潤濕面積w與板的總面積t的比值,如式(10)。

    式中,潤濕面積的是根據(jù)Lan等[32]的研究,以l=0.5處為氣液相界面,利用CFD-post軟件對其進行積分得到。若表示液相在板上的流型為溪流或滴流,表示液相在板上以完全膜狀進行 流動。

    數(shù)為表面張力和慣性力的比值,用數(shù)來表征液體流型的變化可以反映出表面張力對液相分布的影響。從圖4可以看出,表面張力較低的液體(如乙醇、丙酮)在傾斜波紋板上以完全膜狀流進行流動,潤濕面積達到最大。在相同流率下,隨著表面張力的增加,液相在板上的潤濕面積逐漸減小,開始以溪流形式進行流動。另外可以明顯地看到溪流液膜表面呈現(xiàn)一定的波動,在液膜積聚處有明顯的突起,隨后液體又在板上有一定的分散。對于表面張力相當?shù)奈锵?,如甘油和水以不同體積比混合的溶液,隨黏度增加,潤濕面積也逐漸增大。

    3.2 計算潤濕面積新模型

    為了更好地描述表面張力和黏度對液相在傾斜波紋板上分布的影響及其影響程度,在對不同液相物系在規(guī)整填料片上的潤濕面積的計算和分析基礎上,猜想液相的表面張力和黏度對液相在傾斜波紋板上的潤濕面積的影響可以用式(11)進行表述。

    為了確定參數(shù)1、2、3的具體數(shù)值,又對不同表面張力、不同黏度的物系在不同流量在傾斜波紋板上的潤濕情況進行了數(shù)值模擬,最后利用Matlab數(shù)學軟件進行多元非線性最小二乘法擬合,得到各個參數(shù)的具體數(shù)值,如式(12)。

    (12)

    圖5為數(shù)值模擬實驗數(shù)據(jù)點與CFD關聯(lián)式(12)之間的關系圖,計算誤差在8%以內,二者吻合性較好,該模型可以很好地預測液體在傾斜波紋板上的潤濕情況。從作者課題組提出的這個計算潤濕面積新模型中可以看到,液體表面張力和黏度對液相在傾斜波紋板上潤濕面積的影響是相反的,潤濕面積隨表面張力增加而減小,隨黏度的增加而有微量的增加,但表面張力的影響作用要遠大于黏度對潤濕面積的影響。

    3.3 傾斜波紋板上的液膜厚度

    作為液體流動特性重要流體力學參數(shù)之一的液膜厚度,需要對其進行詳細的研究。為了避免液相進口對模擬結果的影響,以距液相進口5mm處的水平面為基準,考察液相物性對液膜厚度的影響。圖6給出了在此水平面上兩組黏度相近但表面張力相差較大的物系之間液膜厚度的比較關系。從圖6中可以看到,低黏度物系組,即低表面張力物系(乙醇)的液膜厚度與高表面張力物系(水)的液膜厚度之間的差值大于高黏度物系組液膜厚度的差值??梢哉f在相同流量下以膜狀流進行流動時,表面張力的增加可以使液膜厚度有微小的增加,但增量不大;在圖6中還可以看到,在壁面附近即圖的左右兩側,受到壁面的影響,液膜厚度急劇增加且遠大于在波紋板面上流動的液膜厚度。因此在工業(yè)生產中要盡量避免壁流現(xiàn)象的出現(xiàn),這會使傳質、傳熱效率在很大程度上降低。

    圖7給出了在同一水平面上兩組表面張力相近但液相黏度相差較大的物系之間液膜厚度的比較關系。從圖7可以看出,兩組物系的液膜厚度均隨液相黏度的增加而增加,但低表面張力物系組液膜厚度增加的幅值小于高表面張力物系組液膜厚度增加的幅值??梢哉f,對高表面張力物系,液相黏度越大,液膜厚度越大。

    為了進一步考察完全膜狀流在沿傾斜波紋板流動過程中液膜厚度的變化情況,給出在完全膜狀流寬度中間位置處的乙醇和丙酮的液相體積分數(shù)云圖以及液膜自由表面、液膜厚度沿波紋板流動方向的變化,如圖8所示。從圖8(a)、8(b)兩個圖中可以看出,液膜自由表面波動趨勢和波紋壁面上的波紋變化幾乎完全重合,即黏度和重力對液膜的流動結構沒有影響,這一數(shù)值模擬結果與谷芳[14][谷芳, 2004 #33]二維模擬的結果一致,可以說液相只有在板面上以完全膜狀流流動時,三維流動模型才可以簡化為二維模型。圖8(c)為膜狀流在沿傾斜波紋板流動過程中液膜厚度的變化情況,受波紋壁面的影響,液膜厚度并不是保持不變的,而是呈現(xiàn)一定的波動性。最大液膜厚度出現(xiàn)在波谷處(即有一定的積液現(xiàn)象),在波峰處液膜厚度最小,且這種不對稱性隨著黏度的增加而增加。

    綜上,根據(jù)傾斜波紋板上液相分布情況的分析研究可知,表面張力和黏度都對液相在傾斜波紋板表面上的分布有影響,表面張力越小,液相在傾斜波紋板片上分布越均勻,潤濕面積越大,液膜寬度越大,液膜厚度越?。浑S黏度的增加,有效相界面積和液膜厚度均增加,但表面張力的影響程度大于黏度的影響程度,為主要影響因素。為此,在研究規(guī)整填料片上液相分布影響因素時,作者重點考察了表面張力對其的影響。數(shù)值模擬過程中取空氣和水為氣液兩相流,僅改變水的表面張力,其他物性保持不變,根據(jù)Elioni等[33]對液相在規(guī)整填料片上接觸角的影響的研究結果,估算了不同表面張力的水在規(guī)整填料片上的接觸角。

    4 表面張力對液相在規(guī)整填料片上分布的影響

    4.1 液相在規(guī)整填料片表面上的分布

    圖9給出了具有不同表面張力的流體在規(guī)整填料片上的流動分布情況,顏色的深淺不僅可以表征流體體積分數(shù)的大小,也可以表征液膜厚度的大小。從圖9中可以看到,表面張力越低,液體在規(guī)整填料表面分布的越均勻;隨著表面張力的逐漸增加,液相在規(guī)整填料片上的流動分布越來越不均勻,溝流現(xiàn)象越來越明顯且對稱性分布也越明顯。另外從圖9中可以很明顯地觀察到,在規(guī)整填料波紋的波峰相交處即填料節(jié)點處有一定的滯液,液相流股在此節(jié)點處混合而后進行再分配,使得液體從波紋板的一個通道流入相鄰通道。為了避免節(jié)點處的積液,在安裝填料的時候可以使相鄰填料片間隔1~2mm,而對于高比表面張力的物系,可以在不改變液相其他物性的前提下加入一定量的表面活性劑,從而使其在填料內的分布更加均勻,以達到強化傳質的效果。

    4.2 液相在規(guī)整填料內的有效面積

    就規(guī)整填料塔而言,有效面積e一般指單位填料體積內液相的相界面積,是影響傳質過程的最重要因素。規(guī)整填料塔內的有效相界面積可以通過物理方法,如電阻法、光傳導以及反射技術測量得到,但在多數(shù)情況下,規(guī)整填料塔內有效相界面積是通過控制快速反應的吸收過程計算測量得到的。Billet等[34]通過整理已有文獻中測量填料塔內潤濕面積的實驗數(shù)據(jù),歸納得出一個適用于求解所有氣液兩相逆流操作的填料塔內的有效相界面積公式,如式(13)。

    Olujic等[35]通過將Onda模型[36]和Delft模型進行合并,得到一個對求解大比表面積,特別是波紋傾角為60°的規(guī)整填料有效相界面積更加有效的公式,如式(14)。

    式中

    (15)

    (17)

    為填料表面上孔的面積占填料比表面積的比值,對Mellapak型填料,。L為液相在規(guī)整填料內的有效流動角度,如式(18)。

    式中,為填料的波紋傾角;為填料兩波峰之間的寬度;為填料波紋峰高。以常見的比表面積為250m2/m3,波紋角度為45°的填料和大氣壓作為參考點,指數(shù)的計算如式(19)。

    (19)

    本文的研究中,CFD數(shù)值模擬求解得到的有效面積與Billet等[34]和Olujic等[35]的求解方式是不相同的。在數(shù)值模擬過程中,利用CFD-post后處理軟件對液相體積分數(shù)進行積分得到有效相界面積。圖10給出了液相在規(guī)整填料內的有效相界面積隨表面張力的變化情況,并將數(shù)值模擬結果與已有經驗公式進行了對比。隨著表面張力的增加,液相在規(guī)整填料內的有效面積逐漸減小。表面張力較低時,CFD數(shù)值模擬結果與Olujic經驗公式[35]較為一致,隨表面張力的增加,二者的差值逐漸增加。這主要是由于一方面Olujic模型[35]中考慮了填料表面穿孔的影響,而該數(shù)值模擬的物理模型表面沒有穿孔。規(guī)整填料表面穿孔不僅有利于液相從一個流道向另一個流道流動,而且有助于形成液滴,增加有效面積。在本次模擬過程中,沒有發(fā)現(xiàn)有液滴的生成。另一方面,Olujic模型[35]的推導是在Onda等[36]散堆填料模型的基礎上利用有機物所做實驗得出的結果,而有機物的表面張力普遍較低,因此Olujic模型更適用于表面張力較低的物系。當液體具有較高表面張力時,CFD數(shù)值模擬結果與Billet模型[34]的預測結果吻合很好,而Billet模型所采用的實驗數(shù)據(jù)中液相的表面張力較高,可以說Billet模型適用于液相表面張力較高的物系。因此,當利用經驗公式來預測規(guī)整填料塔內有效相界面積時,要根據(jù)液相的具體物性選用合適的經驗公式。模擬的結果說明利用CFD數(shù)值模擬不僅可以使液相在規(guī)整填料片上的分布可視化,還能夠準確反應預測液相物性對液相在規(guī)整填料片上分布的影響和有效相界面積值,不受經驗公式適應條件的限制。

    4.3 規(guī)整填料片上液膜厚度

    觀察液相在規(guī)整填料片上的分布,不難發(fā)現(xiàn)液相在規(guī)整填料片上的液膜厚度并不像傾斜波紋板上均勻。在填料節(jié)點處,由于液相的積聚和再分配,使得液膜呈現(xiàn)一定的周期性分布且在節(jié)點附近液膜厚度最大;液膜厚度在徑向和軸向上的分布也都非常不均勻,因此,利用平均液膜厚度來代替局部液膜厚度具有更大的實用價值。研究填料塔的流體力學性能時,雖然可以通過一定的實驗測量出持液量這一重要參數(shù)。但很多情況下是利用有效相界面積和填料內平均液膜厚度進行相應估算得到的。

    多為利用Nusselt液膜理論[11]得到的近似液膜厚度。

    (21)

    式中,p為規(guī)整填料比表面積;t為計算域內液相體積;w為液相在填料片上的潤濕面積;L為液相在規(guī)整填料內按平均液膜厚度的合理性的有效流動角度,用式(18)進行求解得到。

    為了評估用式(22)求解規(guī)整填料平均液膜厚度的合理性,在不同表面張力下,用CFD數(shù)值模擬得到的平均液膜厚度與用Nusselt液膜理論得到的平均液膜厚度進行比較,如圖11所示,該數(shù)值模擬結果與Haroun等[26]的模擬結果相似。

    從圖11中可以看到,表面張力較低時,數(shù)值模擬計算得到的液膜厚度與Nusselt液膜理論得出的液膜厚度值幾乎完全重合,這也說明數(shù)值模擬方法的可行性和準確性;隨表面張力的增加,CFD/Nussel比值也逐漸增加,這說明液體在規(guī)整填料內的平均液膜厚度隨著表面張力的增加而增加。這主要與液相在規(guī)整填料內的分布有關,而液相分布又是表面張力的函數(shù)。表面張力較低時,液相在規(guī)整填料片上分布較均勻,填料內平均液膜厚度較??;表面張力較高時,液相在規(guī)整填料內主要以溝流形式進行流動,平均液膜厚度較高。數(shù)值模擬的結果表明,只有液體在規(guī)整填料內均勻分布時,利用Nusselt液膜理論求解的平均液膜厚度才比較準確,式(21)適用于求解一維平均液膜厚度,液體在規(guī)整填料內以溝流形式進行三維流動時,利用Nusselt液膜理論求解的平均液膜厚度偏小。無論液相在規(guī)整填料內以怎么樣的形式進行流動,利用CFD數(shù)值模擬的方法都可以準確地得到規(guī)整填料內的平均液膜厚度的大小。

    5 結 論

    本文采用計算流體力學中VOF的方法對液相在規(guī)整填料表面上的分布進行了三維數(shù)值模擬。該方法使液相在規(guī)整填料表面的分布可視化,數(shù)值模擬結果與已有文獻中的經驗公式進行了對比,同時給出了液相表面張力和黏度對規(guī)整填料的有效相界面積比和液膜厚度影響程度。

    對于規(guī)整填料的最簡化模型傾斜波紋板,數(shù)值模擬結果表明,液相在傾斜波紋板上的潤濕面積隨表面張力的減小而增加;隨黏度的增加,潤濕面積和液膜厚度都有微量的增加,并提出一個預測板上潤濕面積的新模型。對規(guī)整填料結構單元的數(shù)值模擬表明,來自不同流道的液相流股在填料節(jié)點處發(fā)生混合和再分配現(xiàn)象,規(guī)整填料有效相界面積受表面張力對液相分布的影響而隨表面張力的增加而減??;在填料節(jié)點處,液膜厚度達到最大值,規(guī)整填料各處的液膜厚度值均不相同,無法用Nusselt液膜理論進行準確估算。CFD數(shù)值模擬使液相在規(guī)整填料內的分布可視化,從而可以看到,液相在規(guī)整填料內的流動更傾向于向左右兩邊流動而不是沿著注射點垂直向下流動。

    為了更好地理解規(guī)整填料塔內的基本流程和提高預測模型的準確性,還需要進行更多的研究。作者課題組將進一步利用CFD數(shù)值模擬考察規(guī)整填料的宏觀結構和微觀結構對液相分布和傳質的影響。

    符 號 說 明

    a——波紋板上波紋的波峰,mm ae/ap——有效相界面積比 at——總面積,m2 aw——潤濕面積,m2 B——規(guī)整填料片上峰寬,mm F——體積力源項,m/s2 g——重力加速度,m/s2 h——規(guī)整填料片上峰高,mm n,t——單位法向量和切向量 p——壓力,Pa Re——雷諾數(shù) S——動量源項 u——流體速度,m/s Vt——計算域內液相體積,m3 We——韋伯數(shù) X,Y,Z——坐標,mm α——體積分數(shù) αL——液相在規(guī)整填料內的有效流動角度 β——波紋板的傾角/規(guī)整填料波紋傾角,(°) θ——固液間接觸角,(°) κ——自由表面曲率 λ——波紋板上波紋的波長,mm μ——物系黏度,Pa?s ρ——流體密度,kg/m3 ——Nusselt液膜厚度,mm ——數(shù)值模擬得到的平均液膜厚度,mm σ——表面張力,N/m ω——潤濕面積比率

    [1] Mohamed Ali A,Janses P,Olujic Z. Experimental characterization and CFD simulation of gas distribution performance of liquid redistributors and collectors in packed columns[J].,2003,81:108-115.

    [2] Alekseenko S V,Markovich D M,Evseev A R,et al. Experimental investigation of liquid distribution over structured packing[J].,2008,54:1424-1430.

    [3] Bravo J L,Rocha J A,F(xiàn)air J R. Mass transfer in gauze packings[J].,1985,64:91-95.

    [4] Bravo J L,Rocha J A,F(xiàn)air J R. Pressure drop instructure packings[J].,1986,3:45-49.

    [5] Rocha J A,Bravo J L,F(xiàn)air J R. Distillation columns containing structured packings:A comprehensive model for their performance. 1.Hydraulic models[J].,1993,32:641-651.

    [6] Rocha J A,Bravo J L,F(xiàn)air J R. Distillation columns containing structured packings:A comprehensive model for their performance. 2. Mass transfer model[J].,1996,35:1660-1667.

    [7] Olujic Z. Development of a complete simulation model for predicting the hydraulic and separation performance of distillation columns equipped with structured packings[J].,1997,11:31-46.

    [8] Olujic Z,Kamerbeek A B,de Graauw J. A corrugation geometry based model for efficiency of structured distillation packing[J].,1999,38:683-685.

    [9] 徐崇嗣,樓建中,姜慶泉. 金屬板波填料液流分布的研究[J]. 化工學報,1986,37(4):402-412.

    [10] Dangizer R. Simulating the turbulent mixing process in a static mixing element[J].,1978,60:1-10.

    [11] Nusselt W. Die oberflaechenkondensateion des wasserdampfes[J].,1916,60:569-575.

    [12] Shetty S A,Cerro R L. Estimation of liquid distribution and mass transfer parameters for ordered packings[C]//Institution of Chemical Engineers Symposium Series,1997a,142:853-864.

    [13] Shetty S,Cerro R L. Fundamental liquid flow correlations for the computaion of desigh parameters for ordered packings[J].,1997b,36:771-783.

    [14] 谷芳. 規(guī)整填料局部流動與傳質的計算流體力學研究[D]. 天津:天津大學,2004.

    [15] 陳江波. 高壓下規(guī)整填料塔的計算傳遞和傳質性能[D]. 天津:天津大學,2006.

    [16] Fernandes J,Lisboa P F,Sim?es P C,et al. Application of CFD in the study of supercritical fluid extraction with structured packing:Wetpressure drop calculations[J].,2009,50:61-68.

    [17] Hosseini S H,Shojaee S,Ahmadi G,et al. Computational fluid dynamics studies of dry and wet pressure drops in structured packings[J].,2012,18:1465-1473.

    [18] Adachi T. Velocity and temperature profiles extending over the liquid and gas phases of two-phase flow falling down vertical plates[J].,2013,51:827-832.

    [19] Raynal L,Royon-Lebeaud A. A multi-scale approach for CFD calculations of gas-liquid flow within large size column equipped with structured packing[J].,2007,62:7196-7204.

    [20] Sun B,He L,Liu B T,et al. A new multi-scale model based on CFD and macroscopic calculation for corrugated structured packing column[J].,2013,59:3119-3130.

    [21] Szulczewska B,Zbicinski I,Gorak A. Liquid flow on structured packing:CFD simulation and experimental study[J].,2003,26:580-584.

    [22] Cooke J J,Armstrong L M,Luo K H,et al. Adaptive mesh refinement of gas-liquid flow on an inclined plane[J].,2014,60:297-306.

    [23] Iso Y,Chen X. Flow transition behavior of the wetting flow between the film flow and rivulet flow on an inclined wall[J].,2011,133:091101-091107.

    [24] Iso Y,Huang J,Kato M,et al. Numerical and experimental study on liquid film flows on packing elements in absorbers for post-combustion CO2capture[J].,2013,37:860-868.

    [25] Rafati S A,Hosseini S H,Shojaee S,et al. CFD studies of pressure drop and increasing capacity in Mellapak plus 752.Y structured packing[J].,2011,34:1402-1412.

    [26] Haroun Y,Raynal L,Alix P. Prediction of effective area and liquid hold-up in structured packings by CFD[J].,2014,92:2247-2254.

    [27] Hosseini S H,Shojaee S,Ahmadi G,et al. Prediction of the effective area in structured packings by computational fluid dynamics[J].,2011,50:10833-10842.

    [28] Dai C N,Li Q S,Chen B H. Pressure drop and mass transfer study in structured catalytic packings[J].,2012,98:78-87.

    [29] Dai C N,Li Q S,Chen B H,et al. Simulation of structured catalytic packings in a bubble-point reactor[J].,2013,100:373-383.

    [30] Hirt C W,Nichols B D. Volume of fluid (VOF) method for the dynamics of free boundaries[J]., 1981,39:201-225.

    [31] Brackbill J U,Kothe D B,Zemach C. A continuum method for modelling surface tension[J].,1992,100:335-354.

    [32] Lan H,F(xiàn)riedrich M,Armaly B F,et al. Simulation and measurment of 3-D shear-driven thin liquid film flow in a duct[J].,2008,29:449-459.

    [33] Elioni M A,F(xiàn)air J R. Liquid flow over textured surfaces. 1. Contact angles[J].,1999,38:284-291.

    [34] Billet R,Schultes M. Predicting mass transfer in packed columns[J].,1993,16:1-9.

    [35] Olujic Z,Behrens M,Colli L,et al. Predicting the efficiency of corrugated sheet structured packings with large specific surface area[J].,2004,18:89-96.

    [36] Onda K,Takeuchi H,Okumoto Y. Mass transfer coefficients between gas and liquid phases in packed columns[J].,1968(1):56-62.

    Computational fluid dynamics study of liquid distribution on structured packing surface

    TAN Liyuan1,YUAN Xigang1,KALBASSI Mohammad Ali2

    (1State Key Laboratory of Chemical Engineering,Tianjin University,Tianjin 300072,China;2Air Products Public Limited Company,Walton on Thames,UK)

    A volume of fluid (VOF) model of the computational fluid dynamics (CFD) method was used to simulate the liquid distribution on the structured packing surface. 3D simulation visualized liquid distribution. The results were used to estimate liquid film thickness and effective interfacial area ratio on the structured packing surface quantitatively,and the results were validated with existing correlations in the literature. Different systems with different surface tension values were used in the simulation and the effect of both surface tension and viscosity of liquid phase on liquid distribution were examined. Based on the simulation a new equation to predict the effective interfacial area ratio was proposed. The simulation results showed that the CFD method can be used as an effective tool to provide information on the details of the gas and liquid flows in complex packing geometries.

    computational fluid dynamics (CFD); visualization; liquid film thickness; effective interfacial area ratio; surface tension; viscosity

    TQ 028

    A

    1000–6613(2015)09–3221–12

    10.16085/j.issn.1000-6613.2015.09.003

    2015-02-03;修改稿日期:2015-04-20。

    譚麗媛(1989—),女,碩士研究生,從事化工分離的研究。聯(lián)系人:袁希鋼,教授,從事化工傳質分離及化工系統(tǒng)工程研究。E-mail yuanxg@tju.edu.cn。

    猜你喜歡
    規(guī)整液膜表面張力
    考慮軸彎曲的水潤滑軸承液膜建模方法
    高空高速氣流下平板液膜流動與破裂規(guī)律
    液膜破裂對PCCS降膜的影響*
    300kt/a硫酸系統(tǒng)規(guī)整填料使用情況簡介
    提高日用玻璃陶瓷規(guī)整度和表面光滑度的處理方法
    佛山陶瓷(2016年11期)2016-12-23 08:50:27
    神奇的表面張力
    小布老虎(2016年4期)2016-12-01 05:46:08
    MgO-B2O3-SiO2三元體系熔渣表面張力計算
    上海金屬(2016年2期)2016-11-23 05:34:45
    電梯的建筑化藝術探索
    大觀(2016年9期)2016-11-16 10:31:30
    基于發(fā)音機制的貪婪自適應語音時長規(guī)整算法
    計算機工程(2015年8期)2015-07-03 12:20:34
    CaF2-CaO-Al2O3-MgO-SiO2渣系表面張力計算模型
    上海金屬(2014年3期)2014-12-19 13:09:06
    成人高潮视频无遮挡免费网站| 中文欧美无线码| 欧美精品亚洲一区二区| 色视频在线一区二区三区| 最近中文字幕高清免费大全6| 久久久精品94久久精品| 国产人妻一区二区三区在| 麻豆精品久久久久久蜜桃| 性色av一级| 日韩一区二区三区影片| 一二三四中文在线观看免费高清| 少妇丰满av| 人妻夜夜爽99麻豆av| 你懂的网址亚洲精品在线观看| 在线观看免费视频网站a站| 国产精品久久久久久久久免| 国产淫片久久久久久久久| 亚洲综合精品二区| 丝袜脚勾引网站| 麻豆精品久久久久久蜜桃| 蜜桃亚洲精品一区二区三区| 亚洲成人一二三区av| 在线免费十八禁| 欧美成人a在线观看| 一边亲一边摸免费视频| 精品国产乱码久久久久久小说| 黄片wwwwww| 亚洲综合色惰| av又黄又爽大尺度在线免费看| 亚洲成人中文字幕在线播放| 午夜福利在线在线| 欧美精品一区二区大全| 日本-黄色视频高清免费观看| 日本免费在线观看一区| 欧美精品一区二区大全| 午夜福利在线在线| 亚洲精品视频女| 男女免费视频国产| 三级经典国产精品| 大又大粗又爽又黄少妇毛片口| 小蜜桃在线观看免费完整版高清| 国产伦理片在线播放av一区| 亚洲国产精品999| 性高湖久久久久久久久免费观看| 国产乱人视频| 黄色一级大片看看| 欧美人与善性xxx| 性高湖久久久久久久久免费观看| 欧美三级亚洲精品| xxx大片免费视频| 国内精品宾馆在线| 国产精品麻豆人妻色哟哟久久| 欧美区成人在线视频| www.色视频.com| 97热精品久久久久久| 精品一区二区三区视频在线| av女优亚洲男人天堂| 亚洲国产色片| 国产无遮挡羞羞视频在线观看| 妹子高潮喷水视频| 国产91av在线免费观看| 九九爱精品视频在线观看| 亚洲四区av| 天堂8中文在线网| 午夜老司机福利剧场| 久久久久久久大尺度免费视频| 久久久久久久久久成人| 日本色播在线视频| 免费观看无遮挡的男女| 乱码一卡2卡4卡精品| 黄色日韩在线| 欧美精品国产亚洲| 老司机影院成人| 麻豆精品久久久久久蜜桃| 最近2019中文字幕mv第一页| 成人无遮挡网站| videossex国产| 国产成人免费无遮挡视频| 欧美成人a在线观看| 成年人午夜在线观看视频| 一级av片app| 国产淫语在线视频| 久久女婷五月综合色啪小说| 自拍偷自拍亚洲精品老妇| 最近最新中文字幕免费大全7| 国产成人免费无遮挡视频| 欧美精品亚洲一区二区| 97在线人人人人妻| 精品酒店卫生间| 亚洲av欧美aⅴ国产| 看非洲黑人一级黄片| 免费久久久久久久精品成人欧美视频 | 国产片特级美女逼逼视频| av不卡在线播放| 国产69精品久久久久777片| 国产老妇伦熟女老妇高清| 国产精品人妻久久久影院| 一本色道久久久久久精品综合| 丰满乱子伦码专区| 精品久久国产蜜桃| av福利片在线观看| 男女啪啪激烈高潮av片| 国产精品成人在线| 十分钟在线观看高清视频www | 免费观看av网站的网址| 99热网站在线观看| 深爱激情五月婷婷| 男人和女人高潮做爰伦理| 中文字幕亚洲精品专区| 99国产精品免费福利视频| 人妻夜夜爽99麻豆av| 国产精品福利在线免费观看| 精品人妻偷拍中文字幕| 深夜a级毛片| 观看av在线不卡| 免费久久久久久久精品成人欧美视频 | 亚洲第一区二区三区不卡| av播播在线观看一区| 免费av中文字幕在线| 婷婷色综合www| 中文字幕免费在线视频6| 91精品一卡2卡3卡4卡| 肉色欧美久久久久久久蜜桃| 国产欧美另类精品又又久久亚洲欧美| 欧美成人a在线观看| 国产欧美日韩一区二区三区在线 | 中文字幕精品免费在线观看视频 | 网址你懂的国产日韩在线| 精品国产露脸久久av麻豆| 久久久成人免费电影| av在线播放精品| 亚洲欧美日韩另类电影网站 | av黄色大香蕉| 久久久久久久久久成人| 国产成人一区二区在线| 全区人妻精品视频| 日本免费在线观看一区| 91aial.com中文字幕在线观看| 只有这里有精品99| 亚洲欧美精品专区久久| 久久99热这里只有精品18| 国产精品精品国产色婷婷| 午夜福利网站1000一区二区三区| 永久免费av网站大全| 亚洲一级一片aⅴ在线观看| 日日摸夜夜添夜夜爱| 超碰97精品在线观看| 日韩成人伦理影院| 人体艺术视频欧美日本| 午夜日本视频在线| 男女国产视频网站| 美女高潮的动态| 亚洲在久久综合| 久久国产乱子免费精品| 天堂8中文在线网| 毛片女人毛片| 精品一区二区三卡| 亚洲欧美精品自产自拍| 最新中文字幕久久久久| 国产精品爽爽va在线观看网站| 最近的中文字幕免费完整| 亚洲av在线观看美女高潮| 九九爱精品视频在线观看| 蜜桃久久精品国产亚洲av| 极品教师在线视频| 美女高潮的动态| 高清毛片免费看| av视频免费观看在线观看| 一级毛片黄色毛片免费观看视频| 久久精品久久久久久噜噜老黄| 少妇裸体淫交视频免费看高清| 久久精品国产亚洲网站| av网站免费在线观看视频| 精品视频人人做人人爽| 亚洲丝袜综合中文字幕| 美女高潮的动态| 久久99精品国语久久久| 国产成人精品一,二区| 大陆偷拍与自拍| 人妻夜夜爽99麻豆av| 国内少妇人妻偷人精品xxx网站| 色婷婷av一区二区三区视频| 欧美成人a在线观看| 纯流量卡能插随身wifi吗| 国产精品99久久久久久久久| 日本与韩国留学比较| 久久久精品94久久精品| 中文字幕亚洲精品专区| 男人舔奶头视频| 亚洲,一卡二卡三卡| 亚洲第一av免费看| 一区二区av电影网| 热re99久久精品国产66热6| 大陆偷拍与自拍| 搡老乐熟女国产| 国产美女午夜福利| 熟妇人妻不卡中文字幕| 97在线人人人人妻| a级一级毛片免费在线观看| 永久免费av网站大全| 青春草亚洲视频在线观看| 大香蕉97超碰在线| 精品一区二区三卡| 国产精品国产三级专区第一集| 青青草视频在线视频观看| 国产视频内射| 女人久久www免费人成看片| 高清日韩中文字幕在线| 91精品伊人久久大香线蕉| 最近最新中文字幕大全电影3| 三级国产精品欧美在线观看| 高清av免费在线| 狂野欧美激情性bbbbbb| av专区在线播放| 国产永久视频网站| 国国产精品蜜臀av免费| 国产亚洲午夜精品一区二区久久| 成人18禁高潮啪啪吃奶动态图 | 小蜜桃在线观看免费完整版高清| 国产爱豆传媒在线观看| 九色成人免费人妻av| av不卡在线播放| 99re6热这里在线精品视频| 99久久中文字幕三级久久日本| 久久国内精品自在自线图片| 99国产精品免费福利视频| 大码成人一级视频| 国产在线一区二区三区精| 国产黄色免费在线视频| 男人爽女人下面视频在线观看| 亚洲国产日韩一区二区| 精品久久国产蜜桃| 亚洲国产精品国产精品| 国产黄频视频在线观看| 日韩 亚洲 欧美在线| 成年女人在线观看亚洲视频| 午夜福利在线观看免费完整高清在| 欧美高清成人免费视频www| av在线app专区| 国产在线免费精品| 日本午夜av视频| 久久 成人 亚洲| 久久精品国产a三级三级三级| 最近中文字幕2019免费版| 一本色道久久久久久精品综合| 国产精品熟女久久久久浪| 国产一区二区三区综合在线观看 | 伊人久久国产一区二区| 久久精品国产自在天天线| 狂野欧美激情性bbbbbb| 国产精品一区二区性色av| 91久久精品国产一区二区三区| 18禁动态无遮挡网站| 亚洲欧美精品自产自拍| 少妇猛男粗大的猛烈进出视频| 午夜视频国产福利| 免费不卡的大黄色大毛片视频在线观看| 亚洲国产色片| 黄色日韩在线| 国产av国产精品国产| 大又大粗又爽又黄少妇毛片口| 国产真实伦视频高清在线观看| 一区二区三区乱码不卡18| 中文资源天堂在线| 久久久精品94久久精品| 少妇精品久久久久久久| 日韩不卡一区二区三区视频在线| 嘟嘟电影网在线观看| 制服丝袜香蕉在线| 一区二区三区乱码不卡18| 国产伦理片在线播放av一区| 激情五月婷婷亚洲| 精品一品国产午夜福利视频| av免费在线看不卡| 免费看日本二区| 精品人妻偷拍中文字幕| 亚洲四区av| 国产成人freesex在线| 大片电影免费在线观看免费| 免费黄频网站在线观看国产| 色婷婷av一区二区三区视频| 欧美高清性xxxxhd video| 91午夜精品亚洲一区二区三区| 天天躁日日操中文字幕| 亚洲最大成人中文| 午夜免费鲁丝| 久久亚洲国产成人精品v| www.色视频.com| 国产一级毛片在线| 婷婷色麻豆天堂久久| 五月天丁香电影| 欧美日韩综合久久久久久| 欧美日韩视频高清一区二区三区二| 日韩电影二区| 777米奇影视久久| 中文字幕人妻熟人妻熟丝袜美| 深夜a级毛片| 简卡轻食公司| 欧美xxxx性猛交bbbb| 人人妻人人看人人澡| 黑人高潮一二区| 免费看av在线观看网站| 亚洲精品久久久久久婷婷小说| 久久久久网色| 全区人妻精品视频| 午夜福利影视在线免费观看| 在线天堂最新版资源| 国产成人免费无遮挡视频| 欧美人与善性xxx| 色哟哟·www| 涩涩av久久男人的天堂| 美女国产视频在线观看| 欧美成人午夜免费资源| 看非洲黑人一级黄片| 最近2019中文字幕mv第一页| 1000部很黄的大片| 99久国产av精品国产电影| 性高湖久久久久久久久免费观看| 啦啦啦视频在线资源免费观看| 在线精品无人区一区二区三 | 老司机影院毛片| 国产av国产精品国产| 国产毛片在线视频| 一级毛片久久久久久久久女| 日本欧美视频一区| 亚洲av福利一区| 涩涩av久久男人的天堂| 国产精品国产三级国产av玫瑰| 欧美日韩亚洲高清精品| 哪个播放器可以免费观看大片| 在线免费观看不下载黄p国产| 欧美三级亚洲精品| 一本一本综合久久| 亚洲欧美精品自产自拍| 日本一二三区视频观看| 大话2 男鬼变身卡| 在线 av 中文字幕| 交换朋友夫妻互换小说| 亚洲欧美一区二区三区国产| 亚洲av国产av综合av卡| 亚洲国产色片| 日韩一本色道免费dvd| 亚洲精品国产av成人精品| 国产在线免费精品| 大片免费播放器 马上看| 少妇高潮的动态图| 久久综合国产亚洲精品| a级一级毛片免费在线观看| 国产国拍精品亚洲av在线观看| 女的被弄到高潮叫床怎么办| 日本色播在线视频| 国产一级毛片在线| 国产极品天堂在线| 午夜激情福利司机影院| 日日啪夜夜撸| 免费在线观看成人毛片| 国产精品无大码| 在线免费观看不下载黄p国产| 我的老师免费观看完整版| 国产精品嫩草影院av在线观看| av天堂中文字幕网| 99久久综合免费| 一个人看的www免费观看视频| 欧美zozozo另类| 精品少妇黑人巨大在线播放| 欧美zozozo另类| 国产成人精品福利久久| 王馨瑶露胸无遮挡在线观看| 夜夜看夜夜爽夜夜摸| 亚洲欧美清纯卡通| 国产黄片美女视频| 日韩欧美精品免费久久| 国产av精品麻豆| 一区二区三区乱码不卡18| 插逼视频在线观看| 最后的刺客免费高清国语| 美女脱内裤让男人舔精品视频| 精品久久久精品久久久| 制服丝袜香蕉在线| 狂野欧美激情性bbbbbb| 亚洲欧美精品专区久久| 国产精品久久久久久精品古装| 精品人妻视频免费看| 久久人人爽人人片av| 天天躁夜夜躁狠狠久久av| 亚洲av免费高清在线观看| 寂寞人妻少妇视频99o| 日韩免费高清中文字幕av| 全区人妻精品视频| 亚洲欧美中文字幕日韩二区| 日韩欧美 国产精品| 水蜜桃什么品种好| av在线观看视频网站免费| 99热全是精品| 日本与韩国留学比较| 亚洲精品中文字幕在线视频 | 最近最新中文字幕免费大全7| 亚洲国产精品国产精品| 精品久久久久久电影网| 免费少妇av软件| 国产精品免费大片| 91精品一卡2卡3卡4卡| 久久ye,这里只有精品| 国产毛片在线视频| 卡戴珊不雅视频在线播放| 少妇高潮的动态图| 日韩免费高清中文字幕av| 免费观看在线日韩| 久久久色成人| 这个男人来自地球电影免费观看 | 久久人人爽人人片av| 男女啪啪激烈高潮av片| 久久人人爽人人片av| 午夜福利高清视频| 少妇熟女欧美另类| 国产午夜精品久久久久久一区二区三区| 韩国高清视频一区二区三区| 在线观看一区二区三区| 免费黄频网站在线观看国产| 国产精品欧美亚洲77777| 18+在线观看网站| 午夜免费鲁丝| 成人亚洲欧美一区二区av| 精品久久久久久久末码| 亚洲国产高清在线一区二区三| 麻豆成人午夜福利视频| 日产精品乱码卡一卡2卡三| 国产精品久久久久久精品电影小说 | 久久97久久精品| 国产黄片美女视频| 日韩亚洲欧美综合| 亚洲av国产av综合av卡| 中文字幕亚洲精品专区| 久久久久视频综合| 久久久欧美国产精品| 亚洲av电影在线观看一区二区三区| 久久精品国产亚洲av涩爱| 国产成人a∨麻豆精品| 欧美最新免费一区二区三区| 亚洲精品国产av成人精品| 永久网站在线| 久久久久精品性色| 街头女战士在线观看网站| 亚洲人成网站在线播| 十分钟在线观看高清视频www | 免费看日本二区| 免费黄网站久久成人精品| 老师上课跳d突然被开到最大视频| 婷婷色麻豆天堂久久| 99久国产av精品国产电影| 亚洲色图av天堂| 搡老乐熟女国产| av在线老鸭窝| 男人舔奶头视频| 久久国产亚洲av麻豆专区| 99久久精品一区二区三区| 日本与韩国留学比较| 在现免费观看毛片| 亚洲性久久影院| 我的女老师完整版在线观看| 亚洲精品中文字幕在线视频 | av在线老鸭窝| 日本午夜av视频| 啦啦啦啦在线视频资源| 我的老师免费观看完整版| 日本一二三区视频观看| 欧美区成人在线视频| 18+在线观看网站| 国内揄拍国产精品人妻在线| 妹子高潮喷水视频| 我的老师免费观看完整版| 99久久综合免费| 亚洲成色77777| 国产成人a∨麻豆精品| 国产精品熟女久久久久浪| 午夜福利网站1000一区二区三区| 乱系列少妇在线播放| 22中文网久久字幕| 久久午夜福利片| 网址你懂的国产日韩在线| 欧美xxxx黑人xx丫x性爽| 免费看光身美女| 日韩大片免费观看网站| 欧美成人精品欧美一级黄| 中文天堂在线官网| 少妇猛男粗大的猛烈进出视频| 99久久精品一区二区三区| 久久ye,这里只有精品| 国产伦在线观看视频一区| 亚洲av电影在线观看一区二区三区| 赤兔流量卡办理| 国产精品一及| 中文在线观看免费www的网站| 在线观看三级黄色| 中文字幕久久专区| 丰满人妻一区二区三区视频av| 高清在线视频一区二区三区| 精华霜和精华液先用哪个| 水蜜桃什么品种好| 肉色欧美久久久久久久蜜桃| 午夜福利影视在线免费观看| 欧美最新免费一区二区三区| 深爱激情五月婷婷| 亚洲,欧美,日韩| 亚洲精品国产色婷婷电影| av视频免费观看在线观看| 亚洲国产欧美人成| 国产伦精品一区二区三区四那| 久久人人爽人人爽人人片va| 久久久久久久久久久丰满| 欧美日韩精品成人综合77777| 欧美成人精品欧美一级黄| 亚洲综合色惰| 成人午夜精彩视频在线观看| 国产av国产精品国产| 国产成人精品久久久久久| 少妇 在线观看| 搡女人真爽免费视频火全软件| 亚洲第一av免费看| 国产黄片视频在线免费观看| 肉色欧美久久久久久久蜜桃| 免费看av在线观看网站| 七月丁香在线播放| 亚洲国产精品999| 亚洲va在线va天堂va国产| 国产精品99久久久久久久久| 中文字幕精品免费在线观看视频 | 中文字幕亚洲精品专区| 日本猛色少妇xxxxx猛交久久| 人妻制服诱惑在线中文字幕| 2021少妇久久久久久久久久久| 亚洲美女黄色视频免费看| tube8黄色片| 视频中文字幕在线观看| 三级国产精品欧美在线观看| 亚洲av不卡在线观看| 777米奇影视久久| 日本黄色日本黄色录像| 久久久久网色| 免费观看a级毛片全部| 观看免费一级毛片| 午夜免费鲁丝| 国产女主播在线喷水免费视频网站| 欧美日韩亚洲高清精品| 色哟哟·www| 国产在线视频一区二区| 亚洲精品日本国产第一区| av国产免费在线观看| 中文字幕久久专区| 久久国产亚洲av麻豆专区| 国产乱人视频| 高清在线视频一区二区三区| 亚洲国产精品一区三区| 日本-黄色视频高清免费观看| 久久精品国产亚洲网站| av国产免费在线观看| 久久婷婷青草| 国产免费一区二区三区四区乱码| 汤姆久久久久久久影院中文字幕| 丝袜喷水一区| 免费大片18禁| 国产精品一区二区性色av| 男人狂女人下面高潮的视频| 成人国产av品久久久| 久热久热在线精品观看| 国产男女内射视频| .国产精品久久| 国产黄色视频一区二区在线观看| www.色视频.com| 人妻一区二区av| 亚洲av免费高清在线观看| 国产一区二区在线观看日韩| 青春草视频在线免费观看| 天美传媒精品一区二区| 亚洲精品aⅴ在线观看| 亚洲美女搞黄在线观看| 色哟哟·www| 性高湖久久久久久久久免费观看| 国产中年淑女户外野战色| 色婷婷av一区二区三区视频| 伦精品一区二区三区| 毛片女人毛片| 日韩伦理黄色片| 亚洲人成网站在线观看播放| 精品久久久久久久末码| 黄片无遮挡物在线观看| 国产成人a区在线观看| 久久 成人 亚洲| 2022亚洲国产成人精品| 97在线人人人人妻| 观看美女的网站| 边亲边吃奶的免费视频| 99热这里只有精品一区| 久久影院123| 高清黄色对白视频在线免费看 | 亚洲成色77777| 成年免费大片在线观看| 久久人人爽人人爽人人片va| 国产在线免费精品| 麻豆乱淫一区二区| 最近最新中文字幕免费大全7| 人妻制服诱惑在线中文字幕| 下体分泌物呈黄色| 亚洲国产精品专区欧美| 伦精品一区二区三区| 国产精品熟女久久久久浪| 亚洲欧美成人综合另类久久久| 亚洲国产精品999| 亚洲国产色片| 日日撸夜夜添| 在线精品无人区一区二区三 | 少妇被粗大猛烈的视频| 国产淫语在线视频| 精品久久久精品久久久| 美女主播在线视频|