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

    乙烯氣相聚合流化床反應(yīng)器內(nèi)Geldart B類和Geldart D類顆粒流動(dòng)特性的三維數(shù)值模擬

    2016-06-24 06:49:02車煜田洲張瑞高宇新鄒恩廣王斯晗劉柏平華東理工大學(xué)化學(xué)工程聯(lián)合國家重點(diǎn)實(shí)驗(yàn)室上海007華東理工大學(xué)化工過程先進(jìn)控制與優(yōu)化技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室上海007中國石油石油化工研究院大慶化工研究中心黑龍江大慶674
    化工學(xué)報(bào) 2016年2期
    關(guān)鍵詞:計(jì)算流體力學(xué)流化床聚乙烯

    車煜,田洲,張瑞,高宇新,鄒恩廣,王斯晗,劉柏平(華東理工大學(xué)化學(xué)工程聯(lián)合國家重點(diǎn)實(shí)驗(yàn)室,上海 007;華東理工大學(xué)化工過程先進(jìn)控制與優(yōu)化技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,上海 007;中國石油石油化工研究院大慶化工研究中心,黑龍江 大慶 674)

    乙烯氣相聚合流化床反應(yīng)器內(nèi)Geldart B類和Geldart D類顆粒流動(dòng)特性的三維數(shù)值模擬

    車煜1,田洲2,張瑞3,高宇新3,鄒恩廣3,王斯晗3,劉柏平1
    (1華東理工大學(xué)化學(xué)工程聯(lián)合國家重點(diǎn)實(shí)驗(yàn)室,上海 200237;2華東理工大學(xué)化工過程先進(jìn)控制與優(yōu)化技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,上海 200237;3中國石油石油化工研究院大慶化工研究中心,黑龍江 大慶 163714)

    摘要:乙烯氣相聚合流化床反應(yīng)器的設(shè)計(jì)、操作和優(yōu)化依賴于對(duì)聚合物顆粒粒徑大小和分布、氣泡運(yùn)動(dòng)特性及聚合反應(yīng)狀況的準(zhǔn)確描述。采用Eulerian-Eulerian雙流體模型和群體平衡模型耦合方法對(duì)某乙烯氣相聚合中試規(guī)模的工業(yè)流化床反應(yīng)器分別處于常規(guī)聚合工藝(屬Geldart B類顆粒)和免造粒工藝(屬Geldart D類顆粒)時(shí)床體的氣固流動(dòng)特征以及不同顆粒類型對(duì)反應(yīng)器操作狀態(tài)和顆粒運(yùn)動(dòng)特性的影響進(jìn)行了三維數(shù)值模擬研究。與傳統(tǒng)聚乙烯生產(chǎn)工藝相比,免造粒工藝時(shí)的Geldart D類聚合物顆粒更易聚集于氣體入口處區(qū)域,而且會(huì)產(chǎn)生明顯的旋渦并出現(xiàn)較大的氣泡。研究結(jié)果可為免造粒聚乙烯生產(chǎn)工藝的工業(yè)推廣應(yīng)用提供參考。

    關(guān)鍵詞:流化床;聚乙烯;計(jì)算流體力學(xué);Geldart D類顆粒;免造粒工藝;粒度分布

    2015-07-27收到初稿,2015-09-18收到修改稿。

    聯(lián)系人:田洲,劉柏平。第一作者:車煜(1986—),男,博士研究生。

    Received date: 2015-07-27.

    引 言

    流化床反應(yīng)器因其優(yōu)良的性能和相對(duì)簡(jiǎn)單的結(jié)構(gòu)廣泛應(yīng)用于聚烯烴工業(yè)生產(chǎn)中[1]。在諸多生產(chǎn)工藝中,氣相工藝因具有獨(dú)特的經(jīng)濟(jì)和技術(shù)優(yōu)勢(shì)受到了廣泛關(guān)注,而且隨著技術(shù)的進(jìn)步正在發(fā)揮著越來越重要的作用[2-4]。近年來,節(jié)能降耗已成為化工行業(yè)的迫切需求,乙烯氣相聚合工藝也隨之出現(xiàn)了一些新的技術(shù),如免造粒技術(shù)等[3, 5]。免造粒技術(shù)要求所生產(chǎn)的聚合物顆粒具有特定的粒徑分布,主要包括大粒徑催化劑技術(shù)和工藝過程技術(shù)[6]。目前,工藝過程的研究尚處于起步階段,尤其是針對(duì)大粒徑聚合物顆粒在反應(yīng)器中的流體力學(xué)特性及其對(duì)操作過程的影響等[3]。同時(shí),基于顆粒反應(yīng)器技術(shù)(reactor granule technology,RGT)的高性能聚烯烴材料的開發(fā)也已取得了重要進(jìn)展[7-8],但由于該類聚合物顆粒通常具有較大的粒徑,對(duì)工業(yè)反應(yīng)器的設(shè)計(jì)和操作提出了更高要求。按照流化顆粒粒徑和密度的不同,Geldart將其劃分為4類[9],免造粒工藝中的聚合物顆粒屬于Geldart D類顆粒,傳統(tǒng)聚乙烯生產(chǎn)工藝中的顆粒屬于Geldart B類顆粒。考察在兩種生產(chǎn)工藝下顆粒類型對(duì)于流化床反應(yīng)器操作性能的影響規(guī)律,對(duì)于指導(dǎo)開發(fā)國產(chǎn)免造粒聚乙烯工藝技術(shù)具有重要意義。

    顆粒在反應(yīng)器中的流體動(dòng)力學(xué)行為可以通過計(jì)算流體力學(xué)(CFD)方法進(jìn)行數(shù)值模擬。Khan 等[4]對(duì)生產(chǎn)聚烯烴的流化床反應(yīng)器內(nèi)從氣泡行為特征、氣固相間相互作用和顆粒分布等方面的CFD研究進(jìn)行了綜述。研究還發(fā)現(xiàn)[9-10],不同的顆粒類型對(duì)于反應(yīng)器內(nèi)流場(chǎng)分布具有重要影響,表現(xiàn)為不同的操作狀態(tài)。在聚烯烴工業(yè)生產(chǎn)中,大型流化床反應(yīng)器內(nèi)兩相流體動(dòng)力學(xué)行為對(duì)反應(yīng)器操作的影響更為顯著。同時(shí),為了開發(fā)新的生產(chǎn)工藝或產(chǎn)品牌號(hào)需考察不同類型顆粒對(duì)反應(yīng)器內(nèi)聚合反應(yīng)狀況和氣固流動(dòng)行為的影響。而且,工業(yè)反應(yīng)器的模擬通常也會(huì)產(chǎn)生新的問題(如選擇模型參數(shù)、考察邊界條件、調(diào)整計(jì)算方案等)。目前,工業(yè)規(guī)模烯烴聚合流化床反應(yīng)器的相關(guān)研究較少,尚未見到關(guān)于不同顆粒類型對(duì)烯烴聚合反應(yīng)器操作性能和兩相流動(dòng)特性的公開報(bào)道。

    為了考察傳統(tǒng)聚乙烯生產(chǎn)工藝(Geldart B類顆粒)和免造粒工藝(Geldart D類顆粒)對(duì)流化床反應(yīng)器操作的影響,同時(shí)研究適用于免造粒技術(shù)的操作條件和反應(yīng)器設(shè)計(jì)要求,本工作采用Eulerian-Eulerian雙流體模型和群體平衡模型相耦合的方法,針對(duì)國內(nèi)某中試規(guī)模的氣相法聚乙烯生產(chǎn)工業(yè)流化床反應(yīng)器,研究其處于免造粒工藝操作時(shí)反應(yīng)器內(nèi)兩相流動(dòng)和混合特性以及操作狀態(tài)的變化,同時(shí)與傳統(tǒng)聚乙烯生產(chǎn)工藝進(jìn)行比較,提出能夠滿足免造粒生產(chǎn)技術(shù)要求的操作和設(shè)計(jì)方案。

    1 數(shù)學(xué)模型及其數(shù)值解法

    1.1雙流體模型及顆粒動(dòng)力學(xué)理論

    本研究采用Eulerian-Eulerian雙流體模型描述流化床反應(yīng)器內(nèi)兩相流動(dòng)特征。利用顆粒動(dòng)力學(xué)理論表示固相的應(yīng)力、黏度及壓力等,從而實(shí)現(xiàn)模型方程的封閉。氣固兩相的連續(xù)性方程、動(dòng)量守恒方程、能量守恒方程以及其他與顆粒動(dòng)力學(xué)理論相關(guān)的表達(dá)式參見文獻(xiàn)[4]。

    1.2群體平衡模型及其求解方法

    在多相流動(dòng)和反應(yīng)體系中,常使用群體平衡方程(population balance equation,PBE)描述分散相的大小與分布隨時(shí)間和空間變化的特征[11]。

    在流化床反應(yīng)器中,群體平衡模型可表示為

    式中,x是描述顆??臻g位置的坐標(biāo),u是顆粒群的平均速度,n(L;x,t)是以顆粒粒徑為內(nèi)坐標(biāo)表示的數(shù)密度函數(shù),S(L;x,t)是考慮了顆粒增長(zhǎng)、聚并和破碎行為的源相。

    目前,用來求解PBE的主要方法有矩方法、離散法、加權(quán)殘量法和隨機(jī)法(Monte Carlo法)。其中,矩方法是與CFD耦合求解應(yīng)用最廣泛的方法,分為標(biāo)準(zhǔn)矩方法、積分矩方法和直接積分矩方法等[12]。通常,固體顆粒粒徑分布(particle size distribution,PSD)的kk階矩可以定義為

    式中,m0、m1、m2、m3分別表示單位體積混合物中顆粒的總數(shù)、粒徑、總面積、總體積。

    通過矩變換之后,式(1)可寫為

    與其他矩方法相比,積分矩方法(quadrature method of moments,QMOM)顯示出求解耦合CFD-PBM模型方面的強(qiáng)大優(yōu)勢(shì),廣泛應(yīng)用于多相流研究中[13]。其最早是由McGraw[14]提出,建立在Gaussian積分近似的基礎(chǔ)上,其表達(dá)式為

    其中,加權(quán)(wi)和積分節(jié)點(diǎn)(Li)分別由顆粒PSD的低階矩通過PD(product-difference)算法確定。積分矩方法通常只需求解4~6個(gè)矩便能得到精確的數(shù)值結(jié)果[15]。

    將式(4)代入式(1),并且考慮顆粒增長(zhǎng)時(shí)的PBE可以表示為

    1.3其他模型方程簡(jiǎn)介

    基于前期研究結(jié)果[3],本工作選擇分散相的RNG k-ε湍流模型考察該中試規(guī)模工業(yè)流化床反應(yīng)器內(nèi)的流動(dòng)情況,使用Gidaspow曳力模型描述該反應(yīng)器內(nèi)氣固兩相間的相互作用。為了描述聚合物顆粒的增長(zhǎng)速率,采用Hatzantonis等[16]提出的聚合物顆粒生長(zhǎng)模型。乙烯在Ziegler-Natta催化劑上的聚合反應(yīng)機(jī)理復(fù)雜,包括鏈活化、鏈增長(zhǎng)、鏈轉(zhuǎn)移和鏈?zhǔn)Щ畹确磻?yīng)。為了在CFD-PBM耦合模型中考察乙烯聚合反應(yīng)動(dòng)力學(xué)的影響,結(jié)合Chen等[15]的研究,使用鏈增長(zhǎng)反應(yīng)的動(dòng)力學(xué)模型。

    在乙烯氣相聚合反應(yīng)過程中,氣固相間的傳熱對(duì)于流化床反應(yīng)器的軸、徑向溫度分布和操作狀態(tài)具有重要影響[15],相間傳熱模型可表示為溫度差的函數(shù),如式(6)所示

    其中,hsg=hgs,是兩相傳熱系數(shù)值,可表示為

    其中,Nus為固體顆粒相的Nusselt,可采用應(yīng)用較為廣泛的Ranz-Marshall關(guān)系式表示。

    1.4三維多尺度耦合模型

    本工作采用CFD-PBM耦合模型進(jìn)行該中試規(guī)模工業(yè)流化床反應(yīng)器的數(shù)值模擬,前期研究已對(duì)其進(jìn)行了介紹[3-17],采用Eulerian-Eulerian雙流體模型描述流化床反應(yīng)器中兩相流體動(dòng)力學(xué)行為,而對(duì)于其中具有多分散特性的固體顆粒則使用群體平衡模型進(jìn)行描述。為細(xì)致研究?jī)上嘀g的相互作用和流體力學(xué)行為,采用三維CFD-PBM模型進(jìn)行模擬研究。聚合物顆粒增長(zhǎng)模型以及乙烯聚合反應(yīng)動(dòng)力學(xué)模型等通過用戶自定義函數(shù)(user defined functions,UDF)編譯,然后與CFD-PBM模型進(jìn)行耦合求解。

    首先由CFD模型求解出固含率、固體顆粒速度和溫度等參數(shù),再利用積分矩方法求解PBE中的矩方程,然后利用聚合物顆粒粒徑分布矩求出索特直徑,接著在CFD模型中求出兩相之間的相互作用力,最后對(duì)固含率、顆粒速度和溫度等進(jìn)行更新,再次求解。通過以上步驟,可以實(shí)現(xiàn)三維CFD-PBM模型的耦合求解。本工作采用該模型描述中試規(guī)模的工業(yè)流化床反應(yīng)器,考慮了中試規(guī)模反應(yīng)器建模所帶來的諸如模型參數(shù)選擇、操作條件考察、邊界條件獲取以及求解方案調(diào)整等新問題,研究了兩種不同的聚乙烯生產(chǎn)工藝(常規(guī)工藝的B類和免造粒工藝的D類顆粒體系)下氣固兩相流體動(dòng)力學(xué)特征。

    1.5計(jì)算過程

    該模擬工作是在Fluent 14.0(Ansys Inc., USA)上采用三維模型以雙精度模式進(jìn)行的。在數(shù)值方法上采用有限體積法,空間離散采用一階迎風(fēng)格式。為了獲得更加穩(wěn)定和精確的結(jié)果,瞬態(tài)項(xiàng)采用一階隱式離散格式,壓力速度耦合采用PC-SIMPLE算法并進(jìn)行修正。根據(jù)文獻(xiàn)及前期研究結(jié)果[3, 18-19],選擇數(shù)值迭代的收斂標(biāo)準(zhǔn)為1×10?3,以固定時(shí)間步長(zhǎng)取樣。同時(shí),為了保證計(jì)算結(jié)果的準(zhǔn)確性和計(jì)算過程的穩(wěn)定性,計(jì)算步長(zhǎng)在整個(gè)計(jì)算過程中逐漸從1×10?5s調(diào)整到1×10?4s。所有計(jì)算工作持續(xù)60 s,而且發(fā)現(xiàn)經(jīng)過10 s時(shí)計(jì)算已經(jīng)基本達(dá)到穩(wěn)定狀態(tài),因此選擇10 s到60 s之間進(jìn)行時(shí)間平均計(jì)算。

    模擬對(duì)象為兩種聚乙烯生產(chǎn)工藝:一是需造粒的常規(guī)工藝(traditional PE production process,TPPP),顆粒體系屬Geldart B類顆粒;二是采用免造粒技術(shù)的生產(chǎn)工藝(non-pelletizing PE production process,NPPP),顆粒體系屬Geldart D類顆粒。分別研究了上述兩種工藝條件下流化床反應(yīng)器床層流體力學(xué)行為,考慮了連續(xù)分布的PSD、聚合反應(yīng)以及顆粒增長(zhǎng)等因素。兩種工藝的顆粒粒徑分布如圖1所示。

    圖1 兩種乙烯聚合工藝的初始顆粒粒徑分布Fig.1 Length number density of initial polymer particles intwo ethylene polymerization processes

    2 模擬對(duì)象及參數(shù)設(shè)置

    2.1研究對(duì)象介紹

    本工作針對(duì)某乙烯氣相聚合中試規(guī)模工業(yè)流化床反應(yīng)器進(jìn)行三維數(shù)值模擬。圖2(a)表示該流化床反應(yīng)器的三維計(jì)算區(qū)域,其結(jié)構(gòu)參數(shù)如圖所示??蓪⒃摲磻?yīng)器計(jì)算區(qū)域劃分為3個(gè)部分:密相區(qū)流化段(0~2.0 m)、床層過渡段和床體擴(kuò)大段。對(duì)該計(jì)算區(qū)域采用Gambit 2.4.6(Ansys Inc., USA)在笛卡兒坐標(biāo)系下進(jìn)行三維網(wǎng)格劃分,在氣體入口處和近壁面處進(jìn)行局部網(wǎng)格加密處理,以適應(yīng)流場(chǎng)的變化,三維體網(wǎng)格如圖2(b)所示。同時(shí),針對(duì)B類顆粒已進(jìn)行了網(wǎng)格無關(guān)性的檢驗(yàn)工作(網(wǎng)格數(shù)量分別為99440、202710、486794),發(fā)現(xiàn)當(dāng)網(wǎng)格數(shù)量從202710增加到486794的過程中固含率和固體平均速度的變化均不明顯,但是計(jì)算時(shí)間和成本卻大幅增加。因此,在該模擬研究中采用202710個(gè)六面體網(wǎng)格(入口面:節(jié)點(diǎn)數(shù)1449,四邊形單元1398;出口面:節(jié)點(diǎn)數(shù)1449,四邊形單元1398;垂直段壁面:節(jié)點(diǎn)數(shù)9200,四邊形單元9100;過渡段壁面:節(jié)點(diǎn)數(shù)4000,四邊形單元3900;擴(kuò)大段壁面:節(jié)點(diǎn)數(shù)1600,四邊形單元1500。該三維網(wǎng)格總節(jié)點(diǎn)數(shù)目為211554)進(jìn)行模擬計(jì)算。

    圖2 中試乙烯氣相聚合流化床反應(yīng)器Fig.2 Sketch of pilot-plant FBR

    2.2計(jì)算條件和參數(shù)設(shè)置

    在初始情況下,該流化床反應(yīng)器的固定床層高度為1.137 m,同時(shí)設(shè)定固體顆粒的速度為零,氣體入口處氣速均勻分布。按照反應(yīng)器工業(yè)操作條件,氣相的進(jìn)、出口溫度分別設(shè)定為356.6 K和360.6 K。對(duì)于氣相,壁面條件設(shè)定為無滑移邊界,氣體在壁面處的速度為零,而對(duì)于固相,顆粒沿壁面運(yùn)動(dòng),因此邊界條件設(shè)為自由滑移邊界條件。計(jì)算區(qū)域的頂部出口設(shè)為壓力出口邊界,在氣體入口處設(shè)定為速度入口邊界。氣固兩相物性參數(shù)見表1,其他操作條件和模型參數(shù)的選擇列于表2。

    表1 氣固兩相物性參數(shù)[13, 20]Table 1 Physical properties of gas and solid phases[13, 20]

    表2 模擬中用到的其他操作條件和模型參數(shù)Table 2 Model parameters and computational conditions employed in simulations

    3 結(jié)果與討論

    3.1模擬計(jì)算與實(shí)驗(yàn)結(jié)果的比較

    針對(duì)傳統(tǒng)造粒工藝(B類顆粒體系),通過比較工業(yè)流化床反應(yīng)裝置內(nèi)床層壓降和溫度分布的實(shí)測(cè)數(shù)據(jù)與模擬結(jié)果對(duì)建立的耦合模型進(jìn)行驗(yàn)證和評(píng)價(jià)。圖3為該流化床反應(yīng)器中床層壓降和溫度沿床層高度方向的分布情況。由圖可知,模擬結(jié)果與實(shí)驗(yàn)測(cè)量的床層壓降和溫度分布吻合良好。從圖中還可以看出,床層壓降隨軸向呈線性分布,床層溫度分布除入口處之外比較均勻,有利于流化床反應(yīng)器中聚合反應(yīng)的正常進(jìn)行。從以上分析可知該三維耦合模型可應(yīng)用于中試規(guī)模工業(yè)流化床反應(yīng)器的模擬。

    圖3 中試流化床反應(yīng)器內(nèi)壓降和溫度沿床層軸向的實(shí)驗(yàn)和模擬結(jié)果比較Fig.3 Model validation of pressure drop and temperature distribution along bed height direction

    3.2Geldart B類顆粒流態(tài)化特性分析

    圖4為傳統(tǒng)聚乙烯生產(chǎn)工藝中床體的固含率沿ZOY平面的分布情況。從圖4(a)平均固含率分布云圖可以看出固體顆粒集中于過渡段壁面附近,因?yàn)樵搮^(qū)域能夠降低固體顆粒的運(yùn)動(dòng)速度,同時(shí)將固體顆粒收集,然后固體顆粒沿壁面向下運(yùn)動(dòng),同時(shí)中部的反應(yīng)氣體攜帶著聚合物顆粒向上運(yùn)動(dòng),在兩相接觸過程中進(jìn)行聚合反應(yīng),實(shí)現(xiàn)熱、質(zhì)傳遞。從圖4(a)中固體顆粒的濃度分布可以觀察到其呈現(xiàn)出以固含率表示的典型的環(huán)核流動(dòng)結(jié)構(gòu)[23],關(guān)于該流動(dòng)結(jié)構(gòu)的詳細(xì)描述如圖5所示。圖4(b)是瞬時(shí)固含率分布云圖,可以觀察到氣泡運(yùn)動(dòng)及其在床層表面處破碎以及顆粒沿著壁面回流的現(xiàn)象。

    圖4 Geldart B類顆粒在流化床反應(yīng)器中平均固含率和瞬時(shí)固含率的分布云圖Fig.4 Mean solid holdups and instantaneous solid volume fraction of Geldart B particles in FBR

    圖5(a)為該流化床在不同高度處的固含率分布。由圖可知,不同高度處的固含率數(shù)值除壁面附近外基本一致,說明該區(qū)域聚合物分布比較均勻。而在壁面附近,由于受到上部過渡段向下流動(dòng)顆粒的影響,固含率數(shù)值均有所上升。圖5(b)表示聚合物顆粒在不同床層高度處的平均速度沿徑向的分布情況。由圖可知,床層核區(qū)由于反應(yīng)氣體裹挾著顆粒一起向上運(yùn)動(dòng),速度較大,尤其是沿床層高度方向平均速度呈現(xiàn)出逐漸增加的趨勢(shì),這是因?yàn)闅馀莸木鄄?huì)影響固體顆粒的運(yùn)動(dòng),氣泡向上運(yùn)動(dòng)過程中氣速會(huì)逐漸增大,因此顆粒平均速度也會(huì)逐漸增大。床層壁面處的顆粒也表現(xiàn)出類似的平均速度分布,但是速度值明顯減小,主要是在壁面附近顆粒受到重力作用向下運(yùn)動(dòng)(速度值為負(fù)),如圖5(c)所示。r=0.75 m處,各個(gè)床層高度的顆粒速度值基本保持不變,是顆粒向上運(yùn)動(dòng)(速度值為正)和向下運(yùn)動(dòng)(速度值為負(fù))的分界點(diǎn),因此通常將流化床中兩相流動(dòng)區(qū)域劃分為壁面附近的環(huán)區(qū)和中心區(qū)域的核區(qū)。在整個(gè)流化段,其環(huán)區(qū)和核區(qū)的寬度基本保持不變,兩相流動(dòng)行為比較穩(wěn)定。圖5(c)是顆粒沿軸向的速度分布,顆粒在核區(qū)受到氣體的作用力向上運(yùn)動(dòng),在環(huán)區(qū)受到重力的作用向下運(yùn)動(dòng),在不同的床層高度處呈現(xiàn)出不同的速度分布特征。

    圖5 傳統(tǒng)聚乙烯生產(chǎn)工藝中Geldart B類顆粒在流化床反應(yīng)器內(nèi)沿床層徑向分布情況Fig.5 Solid holdups, solid mean velocity and mean Z velocity of Geldart B particles in FBR

    3.3Geldart D類顆粒流態(tài)化特性分析

    圖6是在免造粒工藝中該反應(yīng)器的固含率分布云圖。在圖6(a)中,固體顆粒明顯聚集于床層壁面附近,與圖4(a)相比在流化段壁面附近固含率數(shù)值更高,而且大量顆粒聚集于過渡段壁面,同時(shí)氣體入口處附近也出現(xiàn)了聚合物顆粒的聚集。這是因?yàn)槊庠炝9に囍械木酆衔镱w粒屬于Geldart D類顆粒,其運(yùn)動(dòng)特性和流化行為與傳統(tǒng)聚乙烯生產(chǎn)工藝中的Geldart B類顆粒不同[24-25]。從圖6(b)中也可觀察到該固含率分布特征,與圖4(b)相比入口處附近出現(xiàn)了顆粒聚集區(qū),而且有較大的氣泡存在。

    圖6 Geldart D類顆粒在流化床反應(yīng)器中平均固含率和瞬時(shí)固含率分布云圖Fig.6 Mean solid holdups and instantaneous solid volume fraction of Geldart D particles in FBR

    圖7為聚合物顆粒沿徑向的固含率、平均速度和軸向平均速度的分布情況。在圖7(a)中,床層高度為0.2 m時(shí)固含率數(shù)值較高,而且與圖5(a)相比在不同床層高度處的固含率數(shù)值均稍高于對(duì)應(yīng)高度處的值,這與氣固兩相間的相互作用力有關(guān)。圖7(b)表示平均顆粒速度在不同高度處沿徑向的分布情況,表現(xiàn)出與圖5(b)相類似的速度分布,區(qū)別僅在于顆粒平均速度的數(shù)值,同時(shí)還表現(xiàn)出不同寬度的環(huán)核流動(dòng)結(jié)構(gòu),這也是由較大的固體顆粒和較高的操作氣速引起的。圖7(c)為軸向顆粒速度分布圖,其大小和方向表現(xiàn)出與圖5(c)明顯不同的趨勢(shì),而且氣體入口處附近顆粒軸向速度會(huì)發(fā)生反向,產(chǎn)生旋渦,這是因?yàn)樵跉怏w入口處產(chǎn)生了大的旋渦,造成氣固兩相運(yùn)動(dòng)速度存在較大差異。

    圖7 免造粒工藝中Geldart D類顆粒在流化床反應(yīng)器密相區(qū)沿床層徑向分布情況Fig.7 Solid holdups, solid mean velocity and mean Z velocity of Geldart D particles in FBR

    3.4Geldart B類和Geldart D類顆粒的速度矢量分布特征

    圖8為兩種聚乙烯生產(chǎn)工藝下固體顆粒在床體密相區(qū)的速度矢量分布圖。由圖8(a)可以看出,B類顆粒在核區(qū)由于受到入口氣體的推動(dòng)作用而向上運(yùn)動(dòng),核區(qū)的分布區(qū)域較寬,而在壁面附近因受到重力等的影響向下運(yùn)動(dòng),分布較窄。與圖8(a)不同的是,圖8(b)中床層入口區(qū)域(0~0.5 m)的顆粒表現(xiàn)出不同的運(yùn)動(dòng)特征(形成旋渦),這是由于固體顆粒粒徑較大和操作氣速較高所致,而在0.5 m以上區(qū)域固體顆粒的運(yùn)動(dòng)趨于穩(wěn)定,表現(xiàn)出與圖8(a)相類似的分布特征,只是受到顆粒粒徑和操作氣速大小影響,環(huán)核流動(dòng)結(jié)構(gòu)的環(huán)區(qū)和核區(qū)厚度呈現(xiàn)出明顯差異。

    圖9為Geldart B類和Geldart D類顆粒瞬時(shí)速度矢量分布圖。從圖中可以明顯觀察到由于氣泡的運(yùn)動(dòng)產(chǎn)生的旋渦的差異,結(jié)合圖5和圖7可以看出,B類顆粒體系中旋渦的產(chǎn)生較為均勻且沒有明顯的顆粒聚集區(qū),而D類顆粒體系中旋渦較大且顆粒的聚集行為明顯,這主要是由D類顆粒所具有的容易產(chǎn)生大氣泡和流化狀態(tài)不穩(wěn)定等動(dòng)力學(xué)特性引起的。

    圖8 Geldart B類和Geldart D類顆粒在流化床反應(yīng)器密相區(qū)的平均顆粒速度矢量分布圖Fig.8 Time-averaged particle velocity vector of Geldart B particles and Geldart D particles in dense region of FBR

    為了定量分析Geldart B類和Geldart D類顆粒體系中氣泡大小的差異,利用圖像分析軟件(Image-Pro Plus 6.0)對(duì)圖9中的氣泡分別進(jìn)行面積和直徑的統(tǒng)計(jì)分析,結(jié)果發(fā)現(xiàn):B類顆粒體系中氣泡的平均面積為3399.85 mm2,平均直徑為92.36 mm;D類顆粒體系中氣泡的平均面積為8413.38 mm2,平均直徑為126.28 mm。以上結(jié)果說明顆粒類型和操作氣速對(duì)氣泡的產(chǎn)生和運(yùn)動(dòng)具有明顯影響,從而影響流化床反應(yīng)器的操作性能。

    圖9 Geldart B類和Geldart D類顆粒在流化床反應(yīng)器密相區(qū)的顆粒瞬時(shí)速度矢量分布圖Fig.9 Instantaneous particle velocity vector of Geldart B particles and Geldart D particles in dense region of FBR

    3.5氣體入口處Geldart B類和Geldart D類顆粒分布

    為便于進(jìn)一步確認(rèn)兩種生產(chǎn)工藝下氣體入口處附近的聚合物顆粒運(yùn)動(dòng)和分布差異,考察了床體入口區(qū)域(0~0.3 m)的固體顆粒特性。圖10表示傳統(tǒng)生產(chǎn)工藝下不同高度處的平均固含率分布,除了入口處固含率數(shù)值較小外,其余各處固含率沿徑向和軸向的分布比較均勻,即在整個(gè)入口附近區(qū)域氣固兩相混合均勻,能夠進(jìn)行正常的操作。但是,免造粒工藝中,如圖11所示,在0.10 m和0.15 m處都可以觀察到明顯的固體顆粒的聚集,這對(duì)于工業(yè)流化床反應(yīng)器的操作極為不利,會(huì)影響聚合反應(yīng)過程的正常進(jìn)行。

    圖10 入口處區(qū)域Geldart B類顆粒分布云圖Fig.10 Contours represented by solid holdups of Geldart B particles in bed inlet region

    圖11 入口處區(qū)域Geldart D類顆粒分布云圖Fig.11 Contours represented by solid holdups of Geldart D particles in bed inlet region

    圖12為兩種生產(chǎn)工藝下氣體入口處附近的固含率、平均速度和平均軸向速度的分布。由圖可知,兩者在固含率數(shù)值、平均速度大小和分布特征以及固體顆粒軸向速度方面存在明顯差異,表明兩種操作工藝下流化床反應(yīng)器入口處附近的固體顆粒表現(xiàn)出完全不同的運(yùn)動(dòng)行為和分布特性。

    圖12 Geldart B類和Geldart D類顆粒在流化床中的固含率、顆粒平均速度和平均軸向顆粒速度沿床層徑向分布情況Fig.12 Solid holdups, solid mean velocity and mean Z velocity of Geldart B particles and Geldart D particles along radial direction in FBR

    3.6改善D類顆粒流態(tài)化的操作和設(shè)計(jì)策略

    為了改善該流化床反應(yīng)器處于免造粒工藝時(shí)的操作性能,需要優(yōu)化反應(yīng)器中氣固兩相的運(yùn)動(dòng)行為?;谝陨戏治?,結(jié)合普通工藝的Geldart B類和免造粒工藝的Geldart D類顆粒的流化性質(zhì)[10],可從兩方面加以考慮:一是增大操作氣速,減少固體顆粒在入口處附近的聚集,加強(qiáng)氣固兩相之間的流動(dòng)和混合[3];二是改進(jìn)分布板結(jié)構(gòu)(如采用直流型分布板),設(shè)計(jì)適應(yīng)Geldart D類顆粒正常流態(tài)化的分布板。通過分布板的氣泡射流改善入口處固體顆粒的運(yùn)動(dòng)和分布[10]。綜合評(píng)價(jià),在工業(yè)生產(chǎn)中增大氣速會(huì)增加能耗,而且還會(huì)增加固體顆粒的夾帶量,同時(shí)由于兩相之間強(qiáng)烈的相互作用,也會(huì)引起床層操作的不穩(wěn)定。而設(shè)計(jì)適應(yīng)免造粒技術(shù)的分布板可以更好地實(shí)現(xiàn)Geldart D類顆粒的正常流化,保證反應(yīng)器在免造粒工藝下的正常操作。

    4 結(jié) 論

    采用三維CFD-PBM耦合模型描述中試規(guī)模乙烯氣相聚合流化床反應(yīng)器,研究了傳統(tǒng)生產(chǎn)工藝(Geldart B類顆粒)和免造粒生產(chǎn)工藝(Geldart D類顆粒)中反應(yīng)器內(nèi)氣固流體動(dòng)力學(xué)特性,得到如下結(jié)論。

    (1)與傳統(tǒng)生產(chǎn)工藝中的Geldart B類顆粒相比,免造粒生產(chǎn)工藝中的Geldart D類顆粒傾向于聚集在氣體入口處附近,造成流化狀態(tài)以及操作的不穩(wěn)定,這主要是由固體顆粒的大?。―類顆粒)以及操作氣速較大引起的。

    (2)對(duì)流化床反應(yīng)器氣體入口處區(qū)域顆粒運(yùn)動(dòng)行為的定性和定量分析發(fā)現(xiàn),Geldart D類顆粒會(huì)產(chǎn)生明顯的旋渦,同時(shí)產(chǎn)生大的氣泡。研究結(jié)果可以為免造粒聚乙烯生產(chǎn)工藝的反應(yīng)器操作和設(shè)計(jì)提供參考。

    符號(hào)說明

    Bag,kk——聚并產(chǎn)生率,s?1

    Bbr,kk——破碎產(chǎn)生率,s?1

    cp,g——?dú)庀啾葻崛?,J·mol?1·K?1

    cp,s——固相比熱容,J·mol?1·K?1

    Dag,kk——聚并損失率,s?1

    Dbr,kk——破碎損失率,s?1

    dp——索特顆粒直徑,m

    ds——顆粒直徑,m

    ΔH ——聚合反應(yīng)熱容,kJ·mol?1

    kk ——矩?cái)?shù)

    L,Li,j——顆粒直徑,m

    p ——壓力,Pa

    ps——顆粒相壓力,Pa

    Qgs——?dú)夤滔嚅g熱傳導(dǎo)速率,W·m?3·s?1

    Rgas——?dú)怏w常數(shù),J·mol?1·K?1

    Res——顆粒Reynolds 數(shù)

    T ——溫度,K

    t——流動(dòng)時(shí)間,s

    ug——?dú)怏w進(jìn)口速度,m·s?1

    αg——?dú)庀囿w積分?jǐn)?shù)

    αs——固相體積分?jǐn)?shù)

    ρg——?dú)怏w密度,kg·m?3

    ρs——固體顆粒密度,kg·m?3

    μg——?dú)怏w黏度,Pa·s

    κg——?dú)怏w熱導(dǎo)率,W·m?1·K?1

    κs——固體熱導(dǎo)率,W·m?1·K?1

    下角標(biāo)

    g——?dú)怏w

    s——固體顆粒

    References

    [1] 張麗霞. Unipol 氣相法聚乙烯工藝技術(shù)進(jìn)展 [J]. 合成樹脂及塑料, 2013, 30(4): 70-74. ZHANG L X. Progress in Unipol gas-phase polyethylene process technology [J]. China Synthetic Resin and Plastics, 2013, 30(4): 70-74.

    [2] XIE T Y, MCAULEY K B, HSU J C, et al. Gas phase ethylene polymerization: production processes, polymer properties, and reactor modeling [J]. Industrial & Engineering Chemistry Research, 1994, 33(3): 449-479.

    [3] CHE Y, TIAN Z, LIU Z, et al. CFD prediction of scale-up effect on the hydrodynamic behaviors of a pilot-plant fluidized bed reactor and preliminary exploration of its application for non-pelletizing polyethylene process [J]. Powder Technology, 2015, 278: 94-110.

    [4] KHAN M J H, HUSSAIN M A, MANSOURPOUR Z, et al. CFD simulation of fluidized bed reactors for polyolefin production — a review [J]. Journal of Industrial and Engineering Chemistry, 2014, 20(6): 3919-3946.

    [5] TERANO M, SUEHIRO K, SAGAE T, et al. Studies in Surface Science and Catalysis[M]. Amsterdam: Elsevier, 2006: 7-12.

    [6] WU L, WANKE S E. Handbook of Transition Metal Polymerization Catalysts[M]. New Jersey: John Wiley & Sons Inc., 2010: 231-259.

    [7] TIAN Z, GU X P, FENG L F, et al. An atmosphere-switching polymerization process: a novel strategy to advanced polyolefin materials [J]. AIChE Journal, 2013, 59(12): 4468-4473.

    [8] GALLI P, VECELLIO G. Technology: driving force behind innovation and growth of polyolefins [J]. Progress in Polymer Science, 2001, 26(8): 1287-1336.

    [9] GELDART D. Types of gas fluidization [J]. Powder Technology, 1973, 7(5): 285-292.

    [10] FAN L S, ZHU C. Principles of Gas-Solid Flows[M]. Cambridge: Cambridge University Press, 2005:371-420.

    [11] 李倩, 程景才, 楊超, 等. 群體平衡方程在攪拌反應(yīng)器模擬中的應(yīng)用 [J]. 化工學(xué)報(bào), 2014, 65(5): 1607-1615. LI Q, CHENG J C, YANG C, et al. Application of population balance equation in numerical simulation of multiphase stirred tanks [J]. CIESC Journal, 2014, 65(5): 1607-1615.

    [12] MARCHISIO D L, FOX R O. Computational Models for Polydisperse Particulate and Multiphase Systems[M]. Cambridge: Cambridge University Press, 2013:47-101.

    [13] YAN W C, LUO Z H, LU Y H, et al. A CFD-PBM- PMLM integratedmodel for the gas-solid flow fields in fluidized bed polymerization reactors [J]. AIChE Journal, 2012, 58(6): 1717-1732.

    [14] MCGRAW R. Description of aerosol dynamics by the quadrature method of moments [J]. Aerosol Science and Technology, 1997, 27(2): 255-265.

    [15] CHEN X Z, LUO Z H, YAN W C, et al. Three- dimensional CFD-PBM coupled model of the temperature fields in fluidized-bed polymerization reactors [J]. AIChE Journal, 2011, 57(12): 3351-3366.

    [16] HATZANTONIS H, GOULAS A, KIPARISSIDES C. A comprehensive model for the prediction of particle-size distribution in catalyzed olefin polymerization fluidized- bed reactors [J]. Chemical Engineering Science, 1998, 53(18): 3251-3267.

    [17] CHE Y, TIAN Z, LIU Z, et al. A CFD-PBM model considering ethylene polymerization for the flow behaviors and particle size distribution of polyethylene in a pilot-plant fluidized bed reactor [J]. Powder Technology, 2015, 286: 107-123.

    [18] SUN J Y, WANG J D, YANG Y R. CFD investigation of particle fluctuation characteristics of bidisperse mixture in a gas-solid fluidized bed [J]. Chemical Engineering Science, 2012, 82: 285-298.

    [19] YAN W C, LI J, LUO Z H. A CFD-PBM coupled model with polymerization kinetics for multizone circulating polymerization reactors [J]. Powder Technology, 2012, 231: 77-87.

    [20] KANELLOPOULOS V, DOMPAZIS G, GUSTAFSSON B, et al. Comprehensive analysis of single-particle growth in heterogeneous olefin polymerization: the random-pore polymeric flow model [J]. Industrial & Engineering Chemistry Research, 2004, 43(17): 5166-5180.

    [21] FAN Rong. Computational fluid dynamics simulation of fluidized bed polymerization reactors[D]. Iowa: Iowa State University, 2006.

    [22] WEI L H, YAN W C, LUO Z H. A preliminary CFD study of the gas-solid flow fields in multizone circulating polymerization reactors [J]. Powder Technology, 2011, 214 (1): 143-154.

    [23] LOHA C, CHATTOPADHYAY H, CHATTERJEE P K. Assessment of drag models in simulating bubbling fluidized bed hydrodynamics [J]. Chemical Engineering Science, 2012, 75: 400-407.

    [24] HAN Y, WANG J J, GU X P, et al. Homogeneous fluidization of Geldart D particles in a gas-solid fluidized bed with a frame impeller [J]. Industrial & Engineering Chemistry Research, 2012, 51(50): 16482-16487.

    [25] 張翠宣, 葉京生, 宋繼田. 噴動(dòng)床研究與進(jìn)展 [J]. 化工進(jìn)展, 2002, 21(9): 630-634. ZHANG C X, YE J S, SONG J T. Development and prospect of modified spouted beds [J]. Chemical Industry and Engineering Progress, 2002, 21(9): 630-634.

    DOI:10.11949/j.issn.0438-1157.20151205

    中圖分類號(hào):TQ 316.3

    文獻(xiàn)標(biāo)志碼:A

    文章編號(hào):0438—1157(2016)02—0519—11

    基金項(xiàng)目:國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃項(xiàng)目(2012CB720502);國家自然科學(xué)基金項(xiàng)目(21406061);上海市自然科學(xué)基金項(xiàng)目(14ZR1410600);引智計(jì)劃項(xiàng)目(B08021);中央高?;究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金項(xiàng)目。

    Corresponding author:TIAN Zhou, tianzhou@ecust.edu.cn; Prof. LIU Boping, boping@ecust.edu.cn supported by the National Basic Research Program of China (2012CB720502), the National Natural Science Foundation of China (21406061), the Shanghai Municipal Natural Science Foundation (14ZR1410600), the Program of Introducing Talents of Discipline to Universities (B08021) and the Fundamental Research Funds for the Central Universities.

    3D numerical simulation of flow characteristics for Geldart B and Geldart D particles in gas phase ethylene polymerization fluidized-bed reactor

    CHE Yu1, TIAN Zhou2, ZHANG Rui3, GAO Yuxin3, ZOU Enguang3, WANG Sihan3, LIU Boping1
    (1State Key Laboratory of Chemical Engineering, East China University of Science and Technology, Shanghai 200237, China;2Key Laboratory of Advanced Control and Optimization for Chemical Processes, Ministry of Education, East China University of Science and Technology, Shanghai 200237, China;3Daqing Petrochemical Research Center, Petrochemical Research Institute of PetroChina, Daqing 163714, Heilongjiang, China)

    Abstract:The status of polyethylene (PE) particle size and distribution, bubble generation and movement, and polymerization reaction in gas phase ethylene polymerization fluidized-bed reactor (FBR) is significant for PE production process, reactor design, optimization and control. Based on 3 dimensional (3D) Eulerian-Eulerian two-fluid model combined with a population balance model (PBM), this work aims to explore the two-phase flow characteristics and the effects of traditional PE production process (Geldart B particles) and non-pelletizing PE production process (Geldart D particles) on the operating behaviors in a pilot-plant FBR. The simulation resultsmatch well with the industrial measured pressure drop and temperature data. It is also found that the polymer particles observably concentrated on the bed inlet region for the effects of Geldart D particles and superfical gas velocity. In addition, the obvious vortexes and large bubbles can be clearly observed in the bed height direction. The results could provide foundation for the extension and application of the non-pelletizing PE production process.

    Key words:fluidized-bed; polyethylene; CFD; Geldart D particles; non-pelletizing PE production process; particle size distribution

    猜你喜歡
    計(jì)算流體力學(xué)流化床聚乙烯
    流化床丙烷脫氫反應(yīng)段的模擬及優(yōu)化
    后吸收法交聯(lián)聚乙烯制備及存儲(chǔ)性研究
    電線電纜(2018年2期)2018-05-19 02:03:43
    《工程流體力學(xué)》教學(xué)方法探討
    基于預(yù)條件技術(shù)的風(fēng)力機(jī)葉片計(jì)算方法研究
    關(guān)于循環(huán)流化床鍋爐集控運(yùn)行研究
    廢棄交聯(lián)聚乙烯回收利用研究進(jìn)展
    中國塑料(2015年9期)2015-10-14 01:12:16
    民用飛機(jī)靜壓源位置誤差修正設(shè)計(jì)研究①
    科技資訊(2015年17期)2015-10-09 21:02:59
    ◆ 塑料管
    單沉浸管流化床內(nèi)離散顆粒數(shù)值模擬
    甘肅:《聚乙烯(PE)再生料》將于9月1日實(shí)施
    精品免费久久久久久久清纯| 亚洲国产中文字幕在线视频| 一区二区三区国产精品乱码| 久久香蕉激情| 88av欧美| 国产在线观看jvid| 亚洲全国av大片| 麻豆久久精品国产亚洲av| 两性午夜刺激爽爽歪歪视频在线观看 | 国产精品 国内视频| 国产黄片美女视频| av有码第一页| 在线观看日韩欧美| 亚洲三区欧美一区| 黄色视频不卡| 村上凉子中文字幕在线| 国产精品一区二区精品视频观看| 成人18禁高潮啪啪吃奶动态图| 97超级碰碰碰精品色视频在线观看| 久久天躁狠狠躁夜夜2o2o| 成人免费观看视频高清| 18禁观看日本| 亚洲国产日韩欧美精品在线观看 | 亚洲精品av麻豆狂野| 美女午夜性视频免费| 午夜成年电影在线免费观看| 免费高清在线观看日韩| 午夜免费成人在线视频| 黄网站色视频无遮挡免费观看| 精品国产美女av久久久久小说| 亚洲专区国产一区二区| 精品久久久久久久久久久久久 | 亚洲欧美激情综合另类| 成人免费观看视频高清| 精品国产美女av久久久久小说| www.自偷自拍.com| 亚洲成人免费电影在线观看| 色老头精品视频在线观看| 看黄色毛片网站| 啦啦啦观看免费观看视频高清| 亚洲男人的天堂狠狠| 国产国语露脸激情在线看| 国产精品综合久久久久久久免费| 国产精品乱码一区二三区的特点| 在线av久久热| 19禁男女啪啪无遮挡网站| 每晚都被弄得嗷嗷叫到高潮| 搞女人的毛片| av片东京热男人的天堂| 一二三四社区在线视频社区8| 日韩欧美一区二区三区在线观看| 看黄色毛片网站| 国产aⅴ精品一区二区三区波| www.熟女人妻精品国产| 高清在线国产一区| √禁漫天堂资源中文www| 性欧美人与动物交配| 97超级碰碰碰精品色视频在线观看| 无人区码免费观看不卡| 国产成人一区二区三区免费视频网站| 免费av毛片视频| 一进一出抽搐动态| 午夜a级毛片| 精品国产国语对白av| 在线看三级毛片| 欧美在线一区亚洲| 亚洲精品色激情综合| 亚洲国产精品sss在线观看| 男人舔女人的私密视频| 黄网站色视频无遮挡免费观看| 色婷婷久久久亚洲欧美| 亚洲国产欧美网| 久久 成人 亚洲| а√天堂www在线а√下载| 国产91精品成人一区二区三区| 欧美日韩黄片免| 亚洲av电影不卡..在线观看| 亚洲一区二区三区不卡视频| 国产91精品成人一区二区三区| 村上凉子中文字幕在线| 在线视频色国产色| 大型黄色视频在线免费观看| www.自偷自拍.com| 国产主播在线观看一区二区| 巨乳人妻的诱惑在线观看| 成在线人永久免费视频| 久久婷婷成人综合色麻豆| 淫妇啪啪啪对白视频| 一级毛片高清免费大全| 成人18禁在线播放| 深夜精品福利| 国产一卡二卡三卡精品| 国产三级黄色录像| 久久久久久国产a免费观看| 国产一区二区激情短视频| 亚洲熟妇熟女久久| 亚洲国产毛片av蜜桃av| 少妇被粗大的猛进出69影院| 久久狼人影院| 中文字幕人妻丝袜一区二区| 国产一卡二卡三卡精品| 97人妻精品一区二区三区麻豆 | 黄色 视频免费看| 国产又黄又爽又无遮挡在线| 欧美成狂野欧美在线观看| 听说在线观看完整版免费高清| 丰满的人妻完整版| 妹子高潮喷水视频| 高潮久久久久久久久久久不卡| 在线av久久热| 亚洲成国产人片在线观看| 国产色视频综合| 久久午夜亚洲精品久久| 在线观看午夜福利视频| 国产av又大| 久久香蕉激情| 精品欧美一区二区三区在线| 色av中文字幕| 成人三级做爰电影| 特大巨黑吊av在线直播 | 一本大道久久a久久精品| 精品国产超薄肉色丝袜足j| 色精品久久人妻99蜜桃| 精品高清国产在线一区| 久久久国产成人精品二区| 最新在线观看一区二区三区| 热99re8久久精品国产| 免费在线观看黄色视频的| 国产欧美日韩一区二区精品| 亚洲人成电影免费在线| 这个男人来自地球电影免费观看| 禁无遮挡网站| 一区二区三区精品91| av超薄肉色丝袜交足视频| a级毛片在线看网站| 久热爱精品视频在线9| 97超级碰碰碰精品色视频在线观看| 久久精品aⅴ一区二区三区四区| 久久99热这里只有精品18| 搡老岳熟女国产| 亚洲国产毛片av蜜桃av| 欧美激情 高清一区二区三区| 1024香蕉在线观看| 少妇被粗大的猛进出69影院| 久久精品aⅴ一区二区三区四区| 男女午夜视频在线观看| 欧美性长视频在线观看| 欧美日韩亚洲国产一区二区在线观看| 国内久久婷婷六月综合欲色啪| 亚洲一区高清亚洲精品| 欧美成人免费av一区二区三区| 亚洲av电影不卡..在线观看| 亚洲人成电影免费在线| 两个人免费观看高清视频| 亚洲激情在线av| 亚洲av电影在线进入| 国产国语露脸激情在线看| 人人妻人人澡人人看| 首页视频小说图片口味搜索| 天堂影院成人在线观看| 精品无人区乱码1区二区| 老司机福利观看| 国产又黄又爽又无遮挡在线| 九色国产91popny在线| 午夜福利一区二区在线看| 久久久久久久久久黄片| 嫩草影视91久久| 午夜福利高清视频| 曰老女人黄片| 日韩大尺度精品在线看网址| 国产亚洲精品综合一区在线观看 | 麻豆av在线久日| 桃红色精品国产亚洲av| 免费高清视频大片| 91成年电影在线观看| 亚洲精品国产精品久久久不卡| 亚洲人成77777在线视频| 亚洲性夜色夜夜综合| 在线观看午夜福利视频| 男人操女人黄网站| 亚洲免费av在线视频| 国产久久久一区二区三区| 亚洲在线自拍视频| av电影中文网址| 亚洲国产高清在线一区二区三 | 亚洲美女黄片视频| 久久精品人妻少妇| 国产伦在线观看视频一区| e午夜精品久久久久久久| 国内毛片毛片毛片毛片毛片| 可以在线观看毛片的网站| 成人三级做爰电影| 在线观看免费午夜福利视频| 国产成人av教育| 在线看三级毛片| 老司机福利观看| 久久人妻av系列| 日日爽夜夜爽网站| 中文字幕精品免费在线观看视频| 亚洲av五月六月丁香网| 亚洲一区二区三区色噜噜| 国产精品亚洲一级av第二区| 国产精品国产高清国产av| 日本一本二区三区精品| 99国产精品一区二区三区| 老司机午夜福利在线观看视频| 一本一本综合久久| 久久国产亚洲av麻豆专区| 9191精品国产免费久久| 不卡一级毛片| 亚洲自拍偷在线| 在线永久观看黄色视频| 19禁男女啪啪无遮挡网站| 嫩草影院精品99| 久久久久久久久久黄片| 国产视频内射| 亚洲av熟女| 国产真实乱freesex| 给我免费播放毛片高清在线观看| 天堂√8在线中文| 国产午夜精品久久久久久| 免费看十八禁软件| 国产精品亚洲av一区麻豆| 久久青草综合色| 黄色成人免费大全| 天天躁狠狠躁夜夜躁狠狠躁| 免费av毛片视频| 成人午夜高清在线视频 | 女人爽到高潮嗷嗷叫在线视频| 动漫黄色视频在线观看| 国产av一区二区精品久久| 欧美激情高清一区二区三区| 国产国语露脸激情在线看| 给我免费播放毛片高清在线观看| 97超级碰碰碰精品色视频在线观看| 在线观看日韩欧美| 亚洲男人的天堂狠狠| 精品久久蜜臀av无| 欧洲精品卡2卡3卡4卡5卡区| 亚洲第一欧美日韩一区二区三区| 欧美黄色淫秽网站| 两个人看的免费小视频| 国产人伦9x9x在线观看| 国产成人av激情在线播放| 可以在线观看毛片的网站| 50天的宝宝边吃奶边哭怎么回事| 老司机深夜福利视频在线观看| 最好的美女福利视频网| 国产私拍福利视频在线观看| 亚洲黑人精品在线| 欧美中文日本在线观看视频| 高潮久久久久久久久久久不卡| 精品国产国语对白av| 午夜福利在线在线| 十八禁网站免费在线| 久久精品91无色码中文字幕| 美女高潮喷水抽搐中文字幕| 日韩 欧美 亚洲 中文字幕| 亚洲中文字幕一区二区三区有码在线看 | 亚洲欧洲精品一区二区精品久久久| 黑人欧美特级aaaaaa片| 窝窝影院91人妻| 成人欧美大片| 国产亚洲精品综合一区在线观看 | 国产成+人综合+亚洲专区| 亚洲精品中文字幕在线视频| 精品久久久久久,| 91麻豆精品激情在线观看国产| 国产精品日韩av在线免费观看| 亚洲av第一区精品v没综合| 久9热在线精品视频| 免费观看精品视频网站| 国产伦一二天堂av在线观看| 久久久久久免费高清国产稀缺| 丝袜在线中文字幕| 久久久国产欧美日韩av| 日本精品一区二区三区蜜桃| 黄频高清免费视频| 亚洲精品久久国产高清桃花| 这个男人来自地球电影免费观看| 免费在线观看黄色视频的| 最新美女视频免费是黄的| 欧美 亚洲 国产 日韩一| 亚洲 国产 在线| 亚洲一区二区三区不卡视频| 久久精品人妻少妇| 成人亚洲精品一区在线观看| 日本精品一区二区三区蜜桃| 国内揄拍国产精品人妻在线 | 美女高潮到喷水免费观看| 国产精品 国内视频| 两性夫妻黄色片| 亚洲一区高清亚洲精品| 在线播放国产精品三级| 看免费av毛片| 免费无遮挡裸体视频| 成人特级黄色片久久久久久久| 国产高清视频在线播放一区| 久久香蕉国产精品| 桃红色精品国产亚洲av| 精品欧美一区二区三区在线| 日韩一卡2卡3卡4卡2021年| 午夜成年电影在线免费观看| 午夜久久久在线观看| 精品高清国产在线一区| 国产精品亚洲美女久久久| 岛国在线观看网站| 怎么达到女性高潮| 午夜视频精品福利| 一区二区日韩欧美中文字幕| 久久天堂一区二区三区四区| 欧美日韩中文字幕国产精品一区二区三区| 18禁美女被吸乳视频| 欧美性猛交黑人性爽| 成人手机av| 久9热在线精品视频| 亚洲国产欧美网| 村上凉子中文字幕在线| 免费看a级黄色片| 日韩欧美一区视频在线观看| 亚洲人成77777在线视频| x7x7x7水蜜桃| 丁香六月欧美| 老熟妇仑乱视频hdxx| 黄色视频不卡| 男人操女人黄网站| xxx96com| 欧美 亚洲 国产 日韩一| 不卡一级毛片| 天天躁夜夜躁狠狠躁躁| 欧美黄色片欧美黄色片| а√天堂www在线а√下载| 男女视频在线观看网站免费 | 免费女性裸体啪啪无遮挡网站| 亚洲精品久久国产高清桃花| 欧美中文日本在线观看视频| 国产1区2区3区精品| 精品国产超薄肉色丝袜足j| 成人三级黄色视频| a级毛片a级免费在线| 香蕉国产在线看| 熟女少妇亚洲综合色aaa.| 国产免费男女视频| 女人高潮潮喷娇喘18禁视频| 亚洲无线在线观看| 无人区码免费观看不卡| 亚洲 欧美一区二区三区| 日韩中文字幕欧美一区二区| 亚洲,欧美精品.| 十八禁人妻一区二区| 国产单亲对白刺激| 亚洲一区二区三区不卡视频| 免费在线观看影片大全网站| 国产av一区二区精品久久| 黄色视频不卡| 在线观看免费日韩欧美大片| 亚洲天堂国产精品一区在线| 91九色精品人成在线观看| 成人国产一区最新在线观看| 国产三级黄色录像| 在线看三级毛片| 国产免费男女视频| 国产亚洲精品第一综合不卡| 叶爱在线成人免费视频播放| 热99re8久久精品国产| 亚洲天堂国产精品一区在线| 99国产精品一区二区三区| 999久久久精品免费观看国产| 天天躁狠狠躁夜夜躁狠狠躁| 免费观看人在逋| av片东京热男人的天堂| 国产成+人综合+亚洲专区| 波多野结衣av一区二区av| 欧美又色又爽又黄视频| 亚洲精品久久成人aⅴ小说| 性色av乱码一区二区三区2| 婷婷丁香在线五月| 老司机在亚洲福利影院| 18禁国产床啪视频网站| 欧美成人一区二区免费高清观看 | 精品不卡国产一区二区三区| 他把我摸到了高潮在线观看| a在线观看视频网站| 最新美女视频免费是黄的| 亚洲成人久久性| 日韩三级视频一区二区三区| 免费在线观看成人毛片| 99久久99久久久精品蜜桃| 国产99白浆流出| 久久中文字幕一级| 国产精品美女特级片免费视频播放器 | 夜夜躁狠狠躁天天躁| 欧美激情久久久久久爽电影| 999久久久精品免费观看国产| 18禁黄网站禁片免费观看直播| 色综合站精品国产| 欧美黑人巨大hd| 国产av一区在线观看免费| 丝袜在线中文字幕| 国产精品二区激情视频| 国产国语露脸激情在线看| 精品久久久久久久人妻蜜臀av| 人人妻人人澡欧美一区二区| 亚洲欧美一区二区三区黑人| 日韩欧美在线二视频| 99久久久亚洲精品蜜臀av| 少妇裸体淫交视频免费看高清 | а√天堂www在线а√下载| 色尼玛亚洲综合影院| 国产亚洲精品综合一区在线观看 | 精品国内亚洲2022精品成人| 丝袜在线中文字幕| 搡老岳熟女国产| 欧美日韩乱码在线| 禁无遮挡网站| 国产成人精品久久二区二区91| 久久久久九九精品影院| 老司机午夜十八禁免费视频| 看免费av毛片| 亚洲真实伦在线观看| av在线天堂中文字幕| 欧美一级毛片孕妇| 国产精品免费视频内射| 亚洲国产欧洲综合997久久, | 国产99白浆流出| 精品久久久久久久毛片微露脸| 91麻豆av在线| 妹子高潮喷水视频| 淫秽高清视频在线观看| 久久久国产成人精品二区| 欧美 亚洲 国产 日韩一| 久久精品aⅴ一区二区三区四区| 黑丝袜美女国产一区| 老司机福利观看| 成人三级黄色视频| 国产视频一区二区在线看| 中文字幕人妻丝袜一区二区| 很黄的视频免费| 国产精品98久久久久久宅男小说| 欧美黑人欧美精品刺激| 欧美日韩福利视频一区二区| 中文资源天堂在线| 亚洲片人在线观看| bbb黄色大片| 国产黄片美女视频| 久久久久久国产a免费观看| 国产激情久久老熟女| 国产精品综合久久久久久久免费| 亚洲熟妇熟女久久| 欧美中文日本在线观看视频| 欧美精品啪啪一区二区三区| 国产又色又爽无遮挡免费看| 中文字幕精品亚洲无线码一区 | 免费无遮挡裸体视频| 久久精品国产亚洲av香蕉五月| 精品久久久久久久末码| 国产又黄又爽又无遮挡在线| 黄网站色视频无遮挡免费观看| 国产精品爽爽va在线观看网站 | 高潮久久久久久久久久久不卡| 美女免费视频网站| 一区二区三区精品91| 看免费av毛片| 精品国产一区二区三区四区第35| 成年人黄色毛片网站| 色婷婷久久久亚洲欧美| 在线免费观看的www视频| 很黄的视频免费| 白带黄色成豆腐渣| 色综合站精品国产| 亚洲第一电影网av| 亚洲精品在线观看二区| 男人的好看免费观看在线视频 | 久久国产亚洲av麻豆专区| 国产精品影院久久| 夜夜爽天天搞| 精品国产乱子伦一区二区三区| 欧美午夜高清在线| 精品久久蜜臀av无| 午夜免费成人在线视频| 91国产中文字幕| 99国产精品一区二区三区| 天堂√8在线中文| 欧美日韩亚洲国产一区二区在线观看| 天天躁夜夜躁狠狠躁躁| 国产成人精品久久二区二区91| 日日夜夜操网爽| 他把我摸到了高潮在线观看| 亚洲一区高清亚洲精品| 日本黄色视频三级网站网址| 国内毛片毛片毛片毛片毛片| 一本精品99久久精品77| 神马国产精品三级电影在线观看 | 免费看美女性在线毛片视频| 男女那种视频在线观看| 伊人久久大香线蕉亚洲五| 中国美女看黄片| 中文字幕精品免费在线观看视频| www.999成人在线观看| 老司机福利观看| 亚洲国产欧美日韩在线播放| 嫩草影院精品99| 日韩欧美国产在线观看| 国产精品久久久av美女十八| 99在线视频只有这里精品首页| 久久久久国产精品人妻aⅴ院| 黄色视频不卡| 久久热在线av| 成人三级黄色视频| 免费在线观看黄色视频的| 一区二区三区激情视频| 亚洲七黄色美女视频| 亚洲 欧美一区二区三区| 国产91精品成人一区二区三区| 狂野欧美激情性xxxx| 久久久久久九九精品二区国产 | 久久精品影院6| 国产精品 欧美亚洲| 欧美+亚洲+日韩+国产| 久久香蕉激情| 老司机在亚洲福利影院| 亚洲性夜色夜夜综合| 日日摸夜夜添夜夜添小说| 1024香蕉在线观看| 欧美成人免费av一区二区三区| 999精品在线视频| 黄频高清免费视频| 国产精品野战在线观看| 国产精品二区激情视频| 法律面前人人平等表现在哪些方面| ponron亚洲| 亚洲av电影在线进入| www.精华液| 色在线成人网| 精品久久久久久久末码| 精品一区二区三区av网在线观看| 一区福利在线观看| 成年女人毛片免费观看观看9| 成人国产一区最新在线观看| 真人一进一出gif抽搐免费| www日本在线高清视频| 亚洲精品色激情综合| 国产激情久久老熟女| 亚洲国产精品999在线| 麻豆成人午夜福利视频| 欧美绝顶高潮抽搐喷水| 亚洲无线在线观看| 国产激情偷乱视频一区二区| 久久中文看片网| 国产精品 欧美亚洲| 在线观看一区二区三区| 黄色片一级片一级黄色片| 精品不卡国产一区二区三区| 亚洲欧美一区二区三区黑人| 真人做人爱边吃奶动态| 久久天堂一区二区三区四区| 美女扒开内裤让男人捅视频| 国产男靠女视频免费网站| 一本一本综合久久| 精品少妇一区二区三区视频日本电影| av在线天堂中文字幕| 日韩大尺度精品在线看网址| 久久中文看片网| 99热这里只有精品一区 | 国产野战对白在线观看| 日本一区二区免费在线视频| 日本成人三级电影网站| 国产精品一区二区三区四区久久 | 日本在线视频免费播放| 亚洲第一av免费看| 两人在一起打扑克的视频| 麻豆成人av在线观看| 亚洲一码二码三码区别大吗| 熟妇人妻久久中文字幕3abv| 国产av不卡久久| 观看免费一级毛片| 国产精品久久久人人做人人爽| 一卡2卡三卡四卡精品乱码亚洲| 国内少妇人妻偷人精品xxx网站 | 婷婷亚洲欧美| 日韩av在线大香蕉| 女人爽到高潮嗷嗷叫在线视频| АⅤ资源中文在线天堂| 国内精品久久久久精免费| 国产精品二区激情视频| 精品久久久久久久久久久久久 | 久久人人精品亚洲av| 日韩有码中文字幕| 午夜视频精品福利| 日本一区二区免费在线视频| 午夜福利高清视频| 亚洲一卡2卡3卡4卡5卡精品中文| 成人av一区二区三区在线看| 国产真实乱freesex| 久久久久久九九精品二区国产 | 国产成人精品无人区| 亚洲午夜精品一区,二区,三区| 欧美zozozo另类| 免费无遮挡裸体视频| 99国产精品一区二区蜜桃av| 色婷婷久久久亚洲欧美| 日日干狠狠操夜夜爽| 天天一区二区日本电影三级| 大型av网站在线播放| 欧美精品亚洲一区二区| 一边摸一边抽搐一进一小说| 欧美成人午夜精品| 99riav亚洲国产免费| 在线观看www视频免费| 精品久久久久久久末码| 少妇 在线观看| a在线观看视频网站| 国产成人欧美| 日本一区二区免费在线视频| 日韩免费av在线播放|