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

    浮頂罐內(nèi)含蠟原油靜態(tài)儲存中的冷卻膠凝規(guī)律

    2017-12-22 05:37:50趙健董航付小明雷啟盟
    化工學(xué)報 2017年12期
    關(guān)鍵詞:含蠟溫降罐底

    趙健,董航,付小明,雷啟盟

    (1東北石油大學(xué)提高采收率教育部重點實驗室,黑龍江 大慶 163318;2大慶油田有限責(zé)任公司儲運銷售分公司,黑龍江 大慶 163000)

    浮頂罐內(nèi)含蠟原油靜態(tài)儲存中的冷卻膠凝規(guī)律

    趙健1,董航1,付小明2,雷啟盟1

    (1東北石油大學(xué)提高采收率教育部重點實驗室,黑龍江 大慶 163318;2大慶油田有限責(zé)任公司儲運銷售分公司,黑龍江 大慶 163000)

    將失流點以下的含蠟原油看作是多孔介質(zhì)體系,以附加比熱容法描述蠟的結(jié)晶潛熱,動量源項方法表征蠟晶網(wǎng)絡(luò)結(jié)構(gòu)對液態(tài)原油的流動阻力,基于有限體積法數(shù)值模擬含蠟原油的冷卻膠凝過程。結(jié)果表明:傳熱機(jī)制和邊界條件主導(dǎo)了凝油結(jié)構(gòu)的演變進(jìn)程。在自然對流作用下,凝油最先在罐底和罐壁所包圍的區(qū)域內(nèi)產(chǎn)生,且其始終是罐內(nèi)膠凝最嚴(yán)重區(qū)域。罐頂最先形成完整的凝油層,其發(fā)展先后經(jīng)歷了慢速增長和快速增長兩個階段,且其凝油層厚度逐漸趨于均勻分布;其次是罐底,其發(fā)展過程與罐頂相反;最后是罐壁,其凝油層的演變具有從罐底沿罐壁向罐頂推進(jìn)的特點。罐內(nèi)對流越強(qiáng),罐頂凝油層的增長速率越緩慢,罐底凝油層的增長速率越快?;跍囟葓黾澳徒Y(jié)構(gòu)的演變規(guī)律,可以將含蠟原油的冷卻過程分為3個階段,即自然對流占主導(dǎo)的第1階段,導(dǎo)熱逐步取代自然對流的第2階段,及以導(dǎo)熱為主導(dǎo)機(jī)制、邊界條件調(diào)控下的第3階段,同時給出了不同階段原油溫度分布和散熱損失規(guī)律的細(xì)節(jié)。

    原油;傳熱;靜態(tài)儲存;膠凝;數(shù)值模擬

    引 言

    目前,我國含蠟原油的總儲存量已達(dá)到近7000萬噸[1],未來隨我國石油戰(zhàn)略儲備工程的逐步完善,其儲存量還將進(jìn)一步增加。與含蠟原油儲存量增大伴隨而來的是業(yè)內(nèi)對其安全、經(jīng)濟(jì)儲存方案的關(guān)注。含蠟原油凝點高、常溫下流動性差,使其在長期儲存中發(fā)生凝罐事故的風(fēng)險升高。準(zhǔn)確掌握含蠟原油罐存過程中的冷卻膠凝規(guī)律,對于制定合理的加熱方案,保障含蠟原油的安全儲存具有重要意義。目前,已有學(xué)者針對相關(guān)問題開展了研究工作,但多數(shù)針對低含蠟原油的儲存?zhèn)鳠徇^程,且在所研究的溫度范圍內(nèi),原油析蠟量很低,主要呈現(xiàn)牛頓流體的流變特征。其中,李旺等[2-4]采用數(shù)值模擬和現(xiàn)場測試方法對原油儲罐內(nèi)的溫度場變化規(guī)律進(jìn)行了研究,探討了太陽輻射、保溫層厚度等對原油傳熱過程的影響規(guī)律。劉佳等[5-6]對原油儲罐內(nèi)的傳熱及流動特性進(jìn)行了研究,考察了儲油高度、加熱蒸氣量、外界環(huán)境溫度及加熱時間對罐內(nèi)原油流場及溫度場的影響。梁文凱等[7]采用FLUENT軟件研究了儲罐中心截面處不同位置的原油溫度,得到了儲罐降溫過程中的溫度分布規(guī)律。陸雅紅等[8]對10×104m3大型立式儲罐的加熱過程進(jìn)行數(shù)值模擬,從熱分析角度確定現(xiàn)有儲罐加熱器是否合理。Cotter等[9-12]以有限體積法數(shù)值求解了原油自然對流溫降過程中的渦量-流函數(shù)形式的控制微分方程組,得到瞬態(tài)自然對流換熱過程中的流動形態(tài)及傳熱特性,并最終給出了一個簡化的溫降計算模型。Vardar[13]數(shù)值分析了一個汽油罐的填充過程,給出了相應(yīng)的流線和溫度等值線圖。Oliveski等[14]采用實驗和數(shù)值模擬方法對模擬罐內(nèi)原油的自然對流換熱過程進(jìn)行了分析。Zhao[15]基于有限元方法對原油儲罐內(nèi)的溫度場分布進(jìn)行了數(shù)值計算。這些研究主要關(guān)注的是邊界條件、儲罐結(jié)構(gòu)和外部環(huán)境等因素對原油傳熱的影響,但對含蠟原油冷卻膠凝規(guī)律的研究相對較少。趙健[16]以附加比熱容方法表征蠟的結(jié)晶潛熱,引入多孔介質(zhì)假設(shè)以動量源項方法表征蠟晶析出對原油流變性的影響,重點分析了原油靜態(tài)儲存過程的溫度場及自然對流演變規(guī)律,初步闡釋了含蠟原油的冷卻膠凝過程。Wang等[1]以內(nèi)熱源表征蠟的結(jié)晶潛熱,采用多孔介質(zhì)的流動阻力表征蠟晶膠凝結(jié)構(gòu)對原油流動的阻礙作用,對膠凝原油經(jīng)加熱轉(zhuǎn)變?yōu)槿苣z過程中的傳熱和流動行為進(jìn)行了深入研究。這些研究雖然實現(xiàn)了對含蠟原油溶膠-凝膠轉(zhuǎn)換過程的模擬,但尚有不少問題,如原油膠凝結(jié)構(gòu)的演變細(xì)節(jié)及其影響因素、膠凝過程對原油傳熱行為的作用規(guī)律等缺少深入的探討。缺失的這些認(rèn)識對剖析含蠟原油長期罐存過程的傳熱特性,闡明其冷卻膠凝規(guī)律具有重要意義。鑒于此,本文以數(shù)值模擬方法對浮頂罐內(nèi)含蠟原油靜態(tài)儲存中的冷卻膠凝過程進(jìn)行研究,獲取含蠟原油由溶膠向凝膠態(tài)轉(zhuǎn)變過程的細(xì)節(jié)特征,探索影響含蠟原油膠凝演變的主導(dǎo)因素及其作用機(jī)制。

    1 問題定義

    目前,浮頂罐已成為大型原油儲庫的首選儲油設(shè)施,因此,以雙盤式浮頂儲油罐為研究對象,為了兼顧計算效率及精度,且更利于突出所研究物理過程的主體特征,特引入以下簡化和假設(shè)條件。

    (1)忽略浮頂上的排水管、支柱及隔板和椽子等支撐結(jié)構(gòu),將其看作由雙層鋼板及所包圍的空氣夾層構(gòu)成,空氣層厚度取為 0.7 m,以當(dāng)量熱導(dǎo)率表征浮頂整體的導(dǎo)熱性能,由于空氣層內(nèi)的對流換熱及部分支撐結(jié)構(gòu)形成的熱橋,浮頂熱導(dǎo)率遠(yuǎn)高于空氣,依據(jù)現(xiàn)場實測數(shù)據(jù)取為1.05 W·m-1·K-1。

    (2)將罐壁看作由等厚度的鋼板和保溫層構(gòu)成,保溫層厚度80 mm,認(rèn)為罐壁的導(dǎo)熱性能一致,并認(rèn)為其與浮頂緊密接觸。由于以儲罐內(nèi)的儲油空間為研究對象,忽略浮頂以上不與原油接觸的罐壁部分。

    (3)將罐底板看作等厚度的鋼板,考慮到罐底部多安裝加熱器,保留罐底部的加熱管。將罐基礎(chǔ)看作由等厚度的瀝青砂、沙土和回填土組成。對于罐底土壤區(qū)域,以恒溫層表征計算模型的軸向邊界,以距離罐壁7 m遠(yuǎn)處為水平絕熱邊界,認(rèn)為該處的土壤溫度不受儲罐影響。

    (4)假定原油儲存過程中儲罐外部環(huán)境條件保持恒定,包括選取恒定的環(huán)境溫度和風(fēng)速,由此更有利于分析特定條件下儲罐內(nèi)凝油結(jié)構(gòu)的演變過程。

    (5)在對罐體結(jié)構(gòu)進(jìn)行簡化的基礎(chǔ)上,忽略原油的周向溫度變化,采用軸對稱模型描述原油的三維傳熱過程。定義儲罐直徑為14.25 m,高度為5 m,儲存過程中環(huán)境溫度保持在-20℃,風(fēng)速7.9 m·s-1,按照上述簡化和假設(shè)條件確定原油罐存過程傳熱系統(tǒng)的計算區(qū)域如圖1所示。

    定義儲罐內(nèi)所儲存原油為含蠟原油,根據(jù)實測數(shù)據(jù)按式(1)擬合其黏溫方程

    將蠟結(jié)晶析出釋放的潛熱包含在原油的比熱容中,以差示掃描量熱法測得含蠟原油的比熱容數(shù)據(jù),并按式(2)擬合其C-T曲線。

    定義原油的熱導(dǎo)率為 0.14 W·m-1·K-1,熱膨脹系數(shù)為 0.000795 K-1,20℃的密度為 870.1 kg·m-3,為了模擬含蠟原油的膠凝過程,定義其滯流點[17]為30℃,失流點為34℃,其他固體介質(zhì)的物性參數(shù)如表1所示。

    表1 固體介質(zhì)的物性參數(shù)Table 1 Physical properties of solid material

    圖1 傳熱系統(tǒng)計算區(qū)域Fig.1 Computation domain of thermal system

    2 數(shù)學(xué)模型及求解算法

    2.1 數(shù)學(xué)模型

    根據(jù)初始油溫和罐體尺寸,預(yù)計溫降初期罐頂內(nèi)壁的 Grashof數(shù)(Gr)可達(dá)到 1×1010,罐壁Gr可達(dá)到1×109,但隨溫降進(jìn)行,原油黏度增大,Gr都將大幅度降低,并且計算中更為關(guān)注的是原油長期儲存的膠凝過程,因此,將罐內(nèi)自然對流看作是層流流動,在文獻(xiàn)[1,18-21]的基礎(chǔ)上,以附加動量源項表征含蠟原油的膠凝行為,建立的數(shù)學(xué)模型如下。

    (1)基本控制方程

    連續(xù)性方程

    軸向動量方程

    徑向動量方程

    其中

    能量方程

    固體導(dǎo)熱微分方程

    式(6)、式(7)中,Su、Sv為動量方程的附加源項,表征因蠟晶結(jié)構(gòu)形成導(dǎo)致的額外流動阻力,這里將失流點以下的含蠟原油看作是多孔介質(zhì)體系,以孔隙度β表征原油的膠凝程度。當(dāng)油溫高于失流點Ts時,認(rèn)為原油完全處于液態(tài),β取值為1,動量源項為0;當(dāng)油溫低于滯流點Tz時,認(rèn)為原油完全膠凝,β取為0,對應(yīng)的動量源項遠(yuǎn)高于自然對流的浮力作用,原油停止流動。油溫從Ts降至Tz時,β=(To-Tz)/(Ts-Tz),按線性規(guī)律由1減小至0,描述原油的逐步膠凝過程。C為常數(shù),用以調(diào)控膠凝過程中流動阻力的演變進(jìn)程。當(dāng)蠟晶網(wǎng)絡(luò)結(jié)構(gòu)對液態(tài)原油的流動阻力超過對流的浮升力時,原油停止流動,此時認(rèn)為含蠟原油的孔隙度為 0,根據(jù)這樣的等量關(guān)系確定C的取值為84280。

    (2)邊界條件

    罐頂外表面

    罐壁外表面

    徑向土壤邊界

    軸向土壤邊界

    計算區(qū)域的對稱軸

    (3)連接條件

    軸向

    徑向

    (4)初始條件

    2.2 數(shù)值求解算法

    基于有限體積法[19,22],在整個計算域上對建立的數(shù)學(xué)模型及邊界條件進(jìn)行離散求解。其中,時間項離散采用一階全隱式離散格式,擴(kuò)散項采用中心差分格式,對流項離散采用具有更高精度的QUICK格式[23]。由于所研究為瞬態(tài)問題,基于 PISO算法[24]實現(xiàn)壓力速度的耦合。所研究的物理現(xiàn)象中,自然對流起主體作用,浮升力影響顯著,采用體積力分?jǐn)?shù)計算處理壓力插值,基于多重網(wǎng)格方法對代數(shù)方程組求解。鑒于計算域具有規(guī)則的幾何結(jié)構(gòu),在進(jìn)行網(wǎng)格剖分時,以結(jié)構(gòu)化四邊形網(wǎng)格離散罐頂、罐壁和罐底基礎(chǔ)及土壤區(qū)域,對臨近原油計算域處進(jìn)行網(wǎng)格加密。對于原油計算域,由于存在加熱管,對計算域進(jìn)行分塊,以非結(jié)構(gòu)化的四邊形網(wǎng)格離散加熱管所在區(qū)域,其余部分以結(jié)構(gòu)化四邊形網(wǎng)格離散,對罐頂、罐壁和罐底內(nèi)表面附近的網(wǎng)格進(jìn)行加密。此過程中,使用的最小網(wǎng)格間距為0.25 mm,以2%的比例增長網(wǎng)格尺寸直到距離計算域邊界 21.15 mm處,此時網(wǎng)格尺寸增大至1.3 mm,而后在軸向和徑向以不同比例逐步增加網(wǎng)格尺寸,延伸至計算區(qū)域的中心位置,最終原油計算域的總網(wǎng)格數(shù)為471459個,網(wǎng)格劃分方案如圖2所示。由于原油儲罐尺寸大,溫降過程緩慢,為了兼顧計算效率及精度,以不同時間步長進(jìn)行瞬態(tài)離散。計算初期,以 0.2 s的步長計算,求解過程穩(wěn)定,且每一步內(nèi)均可快速收斂。因此,逐步提高時間步長至0.5、1、2 s,計算過程均能保持穩(wěn)定,但收斂所需的迭代步數(shù)明顯增多??紤]到隨溫降進(jìn)展,自然對流減弱,計算域內(nèi)的流動及溫度擾動將逐漸減小,計算將更為穩(wěn)定,算法更易收斂,因此,隨計算過程深入,保持以2 s的步長進(jìn)行后續(xù)計算。

    圖2 原油計算域網(wǎng)格劃分方案Fig.2 Grid generation for oil calculation domain

    2.3 驗證過程

    應(yīng)用建立的數(shù)學(xué)模型及求解算法對大慶油田一座5×104m3浮頂儲油罐的溫降過程進(jìn)行了計算,并利用安裝在罐內(nèi)軸向的溫度傳感器對比計算和實測數(shù)據(jù),所選擇的驗證工況中,儲罐初始中心油溫為 39.8℃,液位保持在 11.4 m,平均環(huán)境溫度為-21.2℃,對比結(jié)果如圖3所示。

    由圖3可知,數(shù)值計算結(jié)果展示出了與實測數(shù)據(jù)相同的軸向溫度變化趨勢,且所選取的數(shù)據(jù)中,實測結(jié)果與計算結(jié)果的最大偏差小于0.5℃,相對偏差低于 2%,具有較好的計算精度,兩者偏差的來源主要是由于實測環(huán)境中環(huán)境條件的波動,及計算中選用的原油物性與實際的差異造成。

    圖3 計算結(jié)果與實測數(shù)據(jù)對比Fig.3 Comparison between numerical solutions and test data

    3 結(jié)果分析與討論

    3.1 凝油結(jié)構(gòu)演變過程

    選取初始油溫為40℃,并在該油溫和邊界條件的共同作用下,計算得到初始時刻不同固體介質(zhì)的溫度場分布。在此基礎(chǔ)上,進(jìn)行原油儲存過程的非穩(wěn)態(tài)傳熱計算。選取不同時刻的計算結(jié)果,以溫降過程中含蠟原油的孔隙度分布為依據(jù),對含蠟原油膠凝結(jié)構(gòu)的演變過程進(jìn)行分析。

    圖4展示了儲罐內(nèi)凝油結(jié)構(gòu)的演變進(jìn)程。從圖中可以看出,凝油最先在罐底和罐壁所包圍的區(qū)域內(nèi)產(chǎn)生[圖 4(a)],這說明自然對流在凝油結(jié)構(gòu)的形成和演變過程中具有重要作用,一般而言,罐頂附近的油溫降低最快,但在自然對流作用下,低溫原油要向罐底運移,同樣,罐壁附近的原油降溫后也要沿罐壁向罐底沉積,在兩者的共同作用下,凝油最先在罐底和罐壁所包圍的區(qū)域內(nèi)產(chǎn)生。由圖 4(b)可知,初始凝油結(jié)構(gòu)形成后,其演變過程十分緩慢,55 h內(nèi),其同時沿罐壁和罐底兩個方向推進(jìn),使凝油區(qū)域增大,推動其擴(kuò)展的因素就是來自罐頂和罐壁冷油的沉積作用。隨罐頂油溫降低,原油黏度增大,對流減弱,導(dǎo)致罐頂最先形成能夠覆蓋整個壁面的凝油層。此外,雖然此時罐底沒有凝油層產(chǎn)生,但在自然對流作用下,罐底存在近0.5 m的過渡層,該層原油的孔隙度介于0和1之間,尚沒有完全膠凝,但這為進(jìn)一步凝油層的形成提供了必要的條件。當(dāng)溫降至145 h時,罐底形成了完整的凝油層,且厚度較為均勻,在凝油層上存在過渡層,其范圍與上一時刻相比有明顯增大[圖 4(c)]。同時,罐頂凝油層厚度也明顯增加,但在對流作用下,凝油層厚度分布不均勻。罐壁雖然沒有形成完整的凝油層,但較上一時刻,凝油區(qū)域向罐壁上部擴(kuò)展明顯。溫降至192 h時,罐壁也形成了完整的凝油層,且越靠近罐底,凝油層越厚。儲罐內(nèi)開始形成以凝油層為外殼,中心偏上部區(qū)域為液態(tài)原油,兩者間為過渡層的結(jié)構(gòu)體系。且由圖 4(d)可見,此時,罐底過渡層最厚,凝油層最薄,罐頂和罐壁過渡層很薄,但凝油層較厚,較上一時刻,罐頂不同位置凝油層厚度的差別有所減小。罐內(nèi)凝油結(jié)構(gòu)最強(qiáng)的區(qū)域仍處于罐底和罐壁的結(jié)合部位,這需要引起額外的注意。一般而言,儲罐的進(jìn)出油口位于罐底以上1 m附近的罐壁上,這正接近罐內(nèi)凝油結(jié)構(gòu)最強(qiáng)的區(qū)域,凝油結(jié)構(gòu)的進(jìn)一步擴(kuò)展極有可能阻礙儲罐收發(fā)油過程的正常進(jìn)行。工程上,一般以有限個溫度監(jiān)測點的數(shù)據(jù)作為判斷罐內(nèi)原油儲存安全性的依據(jù),但由于油溫分布不均勻,溫度監(jiān)測點的位置選擇尤為重要,儲罐內(nèi)部進(jìn)出油口附近區(qū)域應(yīng)該作為長期儲存中的重點監(jiān)控位置。隨溫降進(jìn)展,邊界處的凝油結(jié)構(gòu)繼續(xù)向罐中心推進(jìn)[圖 4(e)、(f)],但凝油層的增長速率明顯減小。由于罐內(nèi)對流作用大幅度減弱,罐頂凝油層厚度的差別越來越小,罐壁附近除接近罐頂和罐底的局部區(qū)域外,中心區(qū)域的凝油層厚度已趨于一致。這說明儲罐內(nèi)的傳熱行為逐步由以自然對流換熱占主導(dǎo)轉(zhuǎn)變?yōu)橐詫?dǎo)熱為主導(dǎo)。此外,溫降進(jìn)行到這一階段,罐內(nèi)已經(jīng)沒有完全處于液態(tài)的原油。預(yù)計隨溫降進(jìn)行,罐內(nèi)原油將以導(dǎo)熱為主要傳熱機(jī)制,以邊界條件決定溫度場和凝油結(jié)構(gòu)分布,緩慢進(jìn)入完全的膠凝狀態(tài)。

    圖4 凝油結(jié)構(gòu)演變過程Fig.4 Evolution of gelatinization structure

    3.2 邊界溫度分布及演變過程

    儲罐不同壁面溫度分布的演變對凝油結(jié)構(gòu)的發(fā)展具有重要影響,也體現(xiàn)了罐內(nèi)徑向和軸向油溫分布的演變特征。因此,提取不同時刻罐頂、罐壁和罐底內(nèi)表面的溫度分布數(shù)據(jù)如圖5~圖7所示。由圖可知,罐頂內(nèi)表面溫度與其凝油層厚度的分布規(guī)律相似,低溫位置的凝油層厚度更厚。罐頂內(nèi)表面的溫度分布沿儲罐徑向存在明顯的波動,且越靠近罐壁,這種波動越顯著,充分體現(xiàn)了罐頂自然對流換熱的影響,此時的流場在文獻(xiàn)[12]中有所提及。隨溫降進(jìn)行,罐頂內(nèi)表面溫度同步降低,但中心位置處溫度降低更為明顯,且隨溫度降低,溫度分布的波動幅度逐漸減小,這體現(xiàn)了自然對流的逐步衰退過程。罐壁的溫度分布可以分為3個區(qū)間,靠近罐頂和罐底區(qū)域的溫度明顯更低,中心區(qū)域溫度較高。溫降的前60 h內(nèi),罐壁中心油溫較為均勻,而后隨溫降進(jìn)展,罐壁溫度整體下降,同時靠近罐底的低溫區(qū)域明顯增大,逐步向中心區(qū)域推進(jìn),而靠近罐頂?shù)牡蜏貐^(qū)域擴(kuò)展不明顯,但溫度梯度很大,這與罐壁凝油結(jié)構(gòu)的演變過程相似,也是罐壁自然對流作用的結(jié)果。溫降192 h后,罐壁的中心高溫區(qū)域消失,僅在靠近罐頂處存在一個高溫度點,且隨溫降進(jìn)行,這一高溫點將逐漸消失,罐壁整體溫度分布僅受罐頂和罐底邊界的影響,這一變化過程可以看作是罐壁附近的傳熱行為由自然對流向?qū)徂D(zhuǎn)變的結(jié)果。圖7展示的是罐底的溫度分布,隨溫降進(jìn)行,罐底溫度整體下降,此外,靠近罐壁處形成的低溫區(qū)域也在沿罐底逐漸向罐中心推進(jìn),低溫區(qū)域的覆蓋范圍逐漸增大,這與罐底凝油層厚度的分布規(guī)律一致,且罐底的溫降速率要明顯小于罐頂和罐壁。

    圖5 浮頂內(nèi)表面溫度分布Fig.5 Temperature profile on inner surface of top wall

    圖6 罐壁內(nèi)表面溫度分布Fig.6 Temperature profile on inner surface of sidewall

    圖7 罐底內(nèi)表面溫度分布Fig.7 Temperature profile on inner surface of base wall

    圖8 罐頂凝油層厚度隨時間的變化Fig.8 Thickness of gelatinization layer vs time on top wall

    圖9 罐底凝油層厚度隨時間的變化Fig.9 Thickness of gelatinization layer vs time on base wall

    3.3 凝油層增長速率

    由圖8、圖9可知,罐頂和罐底的凝油層增長過程都可以分為兩個階段,即低速增長和高速增長階段,但其發(fā)展順序不同。對于罐頂,其主體區(qū)域產(chǎn)生凝油結(jié)構(gòu)的時間基本一致,約在41 h。而后,凝油層以較低速度增長,平均速率約為 0.41 mm·h-1。此過程持續(xù)至115 h后,凝油層開始快速增長,其增長速率增加至1.1 mm·h-1。分析其原因,主要是因為溫降初期,罐頂油溫高,自然對流強(qiáng),罐頂附近的低溫原油大量向罐底沉積,延緩了凝油層的增長進(jìn)程。但隨油溫逐漸降低,原油黏度增大,由罐頂至罐底的宏觀對流也隨之減弱,導(dǎo)致低溫原油不再大范圍地向罐底沉積,而是在罐頂內(nèi)表面附近積累,凝油層的增速也明顯加強(qiáng)。與之相反,罐底凝油層的增長過程雖然也可以分為兩個階段,但在主體部分出現(xiàn)凝油后,首先是凝油層的快速增長,而后轉(zhuǎn)變?yōu)槁僭鲩L階段。結(jié)合前述分析可知,罐底凝油層的增長主要是由于罐頂?shù)蜏卦驮谧匀粚α髯饔孟碌某练e導(dǎo)致,但沉積的原油并不會立即形成凝油層,其主要以過渡層的形成存在,這為罐底凝油層的形成提供了基本條件。隨溫降進(jìn)展,這些低溫原油不斷向罐底散熱,油溫降低,并最終在145 h時,開始形成凝油層。凝油層形成初期,其增長速率較快,直至 192 h,其增長速率可達(dá)到 0.46 mm·h-1。但隨溫降進(jìn)行,作為罐底凝油層原料的低溫沉積原油大幅度減小,導(dǎo)致凝油層增長速率降至0.14 mm·h-1。由圖10可知,雖然罐壁不同位置凝油結(jié)構(gòu)的形成時間有所不同,但除接近罐底的部分區(qū)域外,不同位置凝油層的增長速率相差不大,近似為1.4 mm·h-1。根據(jù)圖4可知,罐壁凝油層的演變具有從罐底沿罐壁向罐頂推進(jìn)的特點,這說明其罐壁結(jié)構(gòu)的發(fā)展,除因為罐壁散熱外,自然對流是其主要推動力。

    圖10 罐壁凝油層厚度隨時間的變化Fig.10 Thickness of gelatinization layer vs time on sidewall

    文獻(xiàn)[25-26]采用測試方法對10×104m3浮頂儲油罐的溫度場進(jìn)行測試,依據(jù)安裝在量油孔內(nèi)的18個測溫點的數(shù)據(jù),發(fā)現(xiàn)罐底凝油層增長速率最快,罐頂其次,經(jīng)歷44 d的溫降,罐底凝油層厚度為罐頂?shù)?倍以上,由于量油孔距離罐壁很近,根據(jù)圖4,基本處于儲罐凝油結(jié)構(gòu)最強(qiáng)的區(qū)域,罐底凝油層厚度明顯高于罐頂,但遠(yuǎn)離該區(qū)域,罐底的凝油層厚度要低于罐頂。在其研究中也指明:罐頂凝油層經(jīng)歷了低速增長和快速增長,罐底凝油層經(jīng)歷了快速增長和低速增長兩個階段。這點在數(shù)值模擬的結(jié)果中也得到體現(xiàn),但由于罐體幾何尺寸和油品物性及所處環(huán)境條件不同,不同階段的轉(zhuǎn)換時間,凝油層增長速率還是有所差異。

    3.4 罐體散熱損失變化過程

    由圖11可知,罐頂、罐壁的平均溫度隨儲存時間具有相似的變化規(guī)律,除最初的快速降溫外,整個溫降過程可以分為3個階段,溫度高于30℃時,罐頂和罐壁的溫降曲線近似為兩條平行的直線,基本以相同速率降溫,且不同壁面的溫差也很小,而30℃正是原油的膠凝溫度。此外,該階段的溫降曲線隨時間有明顯的波動,這些都表明溫降的第1階段,自然對流換熱占據(jù)主導(dǎo)作用。當(dāng)壁溫低于30℃后,罐頂、罐壁的降溫速率都明顯增大,這主要是因為壁面膠凝導(dǎo)致熱油的上浮減弱,使壁面溫度持續(xù)下降,此時也是壁面開始有大范圍的凝油層出現(xiàn)的時刻,不同壁面的溫差也開始增加,邊界條件的影響開始顯現(xiàn),這可以看作是溫降的第2階段。此后,隨壁面溫度大幅度下降,壁面與外部環(huán)境的溫差減小,溫降速率又開始降低并最終趨于穩(wěn)定,此為第3階段。相對而言,罐底壁面溫度的變化趨勢明顯不同,溫降初期,在罐內(nèi)宏觀對流作用下,其溫度下降很快,但隨溫降進(jìn)行,對流減弱,罐底的溫降速率又逐漸減小。根據(jù)圖12,不同邊界熱通量也隨儲存時間有顯著變化。由于環(huán)境條件恒定,熱通量主要受邊界油溫和壁面熱阻決定,導(dǎo)致熱通量與壁面溫度具有相同的變化趨勢。但由于罐頂熱阻遠(yuǎn)低于罐壁,導(dǎo)致罐頂熱通量隨儲存時間的變化更為明顯。由于罐底以導(dǎo)熱為主,熱通量最低,也因為土壤的熱容量極大,導(dǎo)致罐底熱通量隨儲存時間的變化很小。也由圖12可知,溫降初期,不同壁面的熱通量差值明顯,罐頂約為罐壁的2倍,考慮到不同邊界的散熱面積,可以確定溫降初期罐頂散熱損失約占總損失的57%,罐壁占39%,罐底占4%。當(dāng)溫降至310 h時,不同壁面熱通量的差值大幅度減小,總散熱損失中罐頂和罐壁各占47%、罐底占6%。

    圖11 邊界平均溫度隨儲存時間的變化Fig.11 Average temperature vs time on boundary

    圖12 邊界平均熱通量隨儲存時間的變化Fig.12 Average heat flow vs time on boundary

    4 結(jié) 論

    (1)傳熱機(jī)制和邊界條件主導(dǎo)了儲罐內(nèi)凝油結(jié)構(gòu)的形成和演變過程。在自然對流作用下,凝油最先在罐底和罐壁所包圍的區(qū)域內(nèi)產(chǎn)生,而后在罐頂形成能夠覆蓋整個壁面的凝油層,其次是罐底,罐壁最后形成凝油層,但罐頂和罐壁的凝油層厚度不均勻。隨溫降進(jìn)行,罐內(nèi)的傳熱機(jī)制逐漸由以自然對流換熱主導(dǎo)轉(zhuǎn)變?yōu)閷?dǎo)熱,邊界條件的影響也越來越明顯,壁面凝油層的厚度也更均勻,儲罐內(nèi)最終形成以凝油層為外殼,中心偏上部區(qū)域為液態(tài)原油,兩者間為過渡層的結(jié)構(gòu)體系。其中,罐底過渡層最厚、凝油層最薄,罐頂和罐壁過渡層最薄、凝油層較厚,但罐內(nèi)凝油結(jié)構(gòu)最強(qiáng)的區(qū)域始終位于罐底和罐壁的結(jié)合部位。

    (2)罐頂凝油層的發(fā)展先后經(jīng)歷了慢速增長和快速增長兩個階段,而罐底凝油層的發(fā)展過程正相反。研究認(rèn)為罐內(nèi)的宏觀對流造成了該種結(jié)果,且對流越強(qiáng),罐頂凝油層的增長速率越慢,罐底凝油層的增長速率越快。除接近罐底的部分區(qū)域外,罐壁不同位置凝油層的增長速率相差不大,且凝油層的演變過程具有從罐底沿罐壁向罐頂推進(jìn)的特點。

    (3)不同壁面溫度分布的演變對凝油結(jié)構(gòu)發(fā)展具有重要影響,也體現(xiàn)了罐內(nèi)徑向和軸向油溫分布的特征。隨溫降進(jìn)行,罐頂內(nèi)表面溫度同步降低,但中心區(qū)域降溫更明顯,溫度分布的波動幅度逐漸減小;罐壁溫度整體下降,但靠近罐底的低溫區(qū)域逐步向中心推進(jìn),而靠近罐頂?shù)牡蜏貐^(qū)域擴(kuò)展不明顯,但溫度梯度很大。罐底除整體存在溫降外,靠近罐壁處形成的低溫區(qū)域沿罐底向罐中心推進(jìn)的特征也十分顯著。

    (4)原油溫降過程可以分為3個階段,即自然對流占主導(dǎo)的第1階段,此階段罐內(nèi)溫降速率基本一致,溫差也很?。浑S溫降進(jìn)展,導(dǎo)熱逐漸取代自然對流主導(dǎo)罐內(nèi)傳熱,凝油結(jié)構(gòu)也開始發(fā)展,罐內(nèi)溫差及邊界溫降速率都明顯增大;而后溫降進(jìn)入第3階段,溫度場在以導(dǎo)熱主導(dǎo),邊界條件調(diào)控下重新分布。在這一過程中,儲罐散熱損失受邊界溫度影響顯著,且隨溫降進(jìn)行,不同邊界熱通量的差異逐漸減小。

    符 號 說 明

    cpo——原油比熱容,J·kg-1·℃-1

    g——重力加速度,m·s-2

    Hjc——罐基礎(chǔ)總厚度,m

    Hsoil——土壤恒溫層深度,m

    Htank——儲罐總高度,m

    hroof——罐頂外表面?zhèn)鳠嵯禂?shù),W·m-2·℃-1

    hsidewall——罐壁外表面?zhèn)鳠嵯禂?shù),W·m-2·℃-1

    p——壓力,Pa

    Rsoil——土壤區(qū)域總外徑,m

    Rtank——儲罐總外徑,m

    r——徑向坐標(biāo)

    Tair——大氣溫度,℃

    Th——土壤恒溫層溫度,℃

    Ti——固體介質(zhì)溫度,℃

    T

    o——原油溫度,℃

    Troof——罐頂外表面溫度,℃

    Tsidewall——罐壁外表面溫度,℃

    Tsoil,r——土壤熱影響區(qū)邊界溫度,℃

    Tsoil,z——土壤恒溫層溫度,℃

    T0——初始油溫,℃

    t——溫降時間,s

    u——軸向速度,m·s-1

    u0——初始軸向速度,m·s-1

    v——徑向速度,m·s-1

    v0——初始徑向速度,m·s-1

    x——軸向坐標(biāo)

    λi——固體介質(zhì)的熱導(dǎo)率,W·m-1·℃-1

    λo——原油熱導(dǎo)率,W·m-1·℃-1

    λroof——罐頂熱導(dǎo)率,W·m-1·℃-1

    λsidewal——保溫層熱導(dǎo)率,W·m-1·℃-1

    λsoil——土壤熱導(dǎo)率,W·m-1·℃-1

    μo——原油動力黏度,Pa·s

    ρo——原油密度,kg·m-3

    下角標(biāo)

    air——大氣環(huán)境

    h——恒溫層

    i——固體介質(zhì)

    jc——罐基礎(chǔ)

    o——原油

    roof——罐頂

    sidewall——罐壁

    soil——土壤

    0——初始時刻

    [1]WANG M,YU G J,ZHANG X Y,et al.Numerical investigation of melting of waxy crude oil in an oil tank[J].Applied Thermal Engineering,2017,115:81-90.

    [2]李旺,王情愿,李瑞龍,等.大型浮頂油罐溫度場數(shù)值模擬[J].化工學(xué)報,2011,62(S1):108-112.LI W,WANG Q Y,LI R L,et al.Numerical study on temperature field of a large floating roof oil tank[J].CIESC Journal,2011,62(S1):108-112.

    [3]李旺.大型浮頂油罐溫度場數(shù)值模擬方法及規(guī)律研究[D].北京:中國石油大學(xué)(北京),2013.LI W.Study on numerical simulation methods and regularities for temperature fi eld of large fl oating roof oil tank [D].Beijing:China University of Petroleum,2013.

    [4]李超,劉人瑋,李旺,等.大型浮頂儲罐原油溫度場實驗測試研究[J].工程熱物理學(xué)報,2013,34(12):2332-2334.LI C,LIU R W,LI W,et al.Experimental study on temperature field of a large scaled floating roof oil tank[J].Journal of Engineering Thermophysics,2013,34(12):2332-2334.

    [5]劉佳,侯磊,陳雪嬌.10×104m3浮頂罐罐壁附近原油溫度分布數(shù)值模擬[J].油氣儲運,2015,34(3):248-253.LIU J,HOU L,CHEN X J.Numerical simulation for the temperature profile near the sidewall of a 10×104m3floating roof oil tank[J].Oil& Gas Storage and Transportation,2015,34(3):248-253.

    [6]劉佳.原油儲罐溫度場數(shù)值模擬及蒸汽加熱節(jié)能降耗研究[D].北京:中國石油大學(xué),2014.LIU J.Numerical simulation for the temperature profile of an oil tank and study on the energy conservation for steam heating [D].Beijing:China University of Petroleum,2014.

    [7]梁文凱,鄧文俊,叢潤芝,等.基于 FLUENT的儲罐內(nèi)原油溫度分布規(guī)律研究[J].遼寧石油化工大學(xué)學(xué)報,2014,34(5):39-43.LIANG W K,DENG W J,CONG R Z,et al.Research on the temperature profile rule of oil tank based on FLUENT[J].Journal of Liaoning University of Petroleum & Chemical Technology,2014,34(5):39-43.

    [8]陸雅紅,吳江濤.原油儲罐盤管式蒸汽加熱器優(yōu)化設(shè)計[J].化學(xué)工程,2010,38(10):69-72.LU Y H,WU J T.Optimization design for the steam heating tube of oil tank[J].Chemical Engineering(China),2010,38(10):69-72.

    [9]COTTER M,CHARLES M.Transient cooling of petroleum by natural convection in cylindrical storage tanks(Ⅰ):Development and testing the numerical simulator[J].International Journal of Heat and Mass Transfer,1993,36:2165-2174.

    [10]COTTER M,CHARLES M.Transient cooling of petroleum by natural convection in cylindrical storage tanks(Ⅱ):Effect of heat transfer coefficient,aspect ratio and temperature-dependent viscosity[J].International Journal of Heat and Mass Transfer,1993,36:2175-2182.

    [11]COTTER M,CHARLES M.Transient cooling of petroleum by natural convection in cylindrical storage tanks:a simplified heat loss model[J].The Canadian Journal of Chemical Engineering,1992,70:1090-1093.

    [12]COTTER M.Transient natural convection in petroleum storage tanks[D].Toronto,Canada:University of Toronto,1991.

    [13]VARDAR N.Numerical analysis of the transient turbulent flow in a fuel oil storage tank[J].International Journal of Heat and Mass Transfer,2003,46:3429-3440.

    [14]OLIVESKI R D C,MACAGNAN M H,COPETTI J B,et al.Natural convection in a tank of oil:experimental validation of a numerical code with prescribed boundary condition[J].Exp.Therm.Fluid Sci.,2005,29:671-680.

    [15]ZHAO B.Numerical simulation for the temperature changing rule of the crude oil in a storage tank based on the wavelet finite element method[J].J.Therm.Anal.Calorim.,2012,107:387-393.

    [16]趙健.高寒地區(qū)原油儲存過程中的傳熱問題研究及工藝方案優(yōu)化[D].大慶:東北石油大學(xué),2014.ZHAO J.Heat-transfer research and process scheme optimization of oil storage in alpine cold region[D].Daqing:Northeast Petroleum University,2014.

    [17]李才,劉云龍,蘇仲勛.管內(nèi)含蠟原油降溫過程中的放熱問題[J].油田地面工程,1994,13(1):18-20.LI C,LIU Y L,SU Z X.The heat transfer problem during the cooling of waxy crude oil in the pipe[J].Oil-Gas Field Surface Engineering,1994,13(1):18-20.

    [18]HE Y L,TAO W Q,QU Z G,et al.Steady natural convection in a vertical cylindrical envelope with adiabatic lateral wall[J].International Journal of Heat and Mass Transfer,2004,47:3131-3144.

    [19]陶文銓.數(shù)值傳熱學(xué)[M].2版.西安:西安交通大學(xué)出版社,2001.TAO W Q.Numerical Heat Transfer[M].2nd ed.Xi’an:Xi’an Jiaotong University Press,2001.

    [20]宇波.流動與傳熱數(shù)值計算——若干問題的研究與探討[M].北京:科學(xué)出版社,2015.YU B.Numerical Simulation of Flow and Heat Transfer—the Research and Discussion for Several Problem[M].Beijing:Science Press,2015.

    [21]LI W,YU B,WANG Y,et al.Study on general governing equations of computational heat transfer and fluid flow[J].Communications in Computational Physics,2012,12(5):1482-1494.

    [22]PATANKAR V.Numerical Heat Transfer and Fluid Flow[M].United States:McGraw Hill Press,1980.

    [23]LEONARD B P,MOKHTARI S.Ultra-sharp nonoscillatory convection schemes for high-speed steady multidimensional flow[R].Cleveland,OH,United States:NASA Lewis Research Center,1990.

    [24]OLIVEIRA P J R,ISSA R I.An improved PISO algorithm for the computation of buoyant driven flows[J].Numerical Heat Transfer,Part B:Fundamentals,2001,40:473-493.

    [25]于達(dá),方徐應(yīng),李東風(fēng),等.大型浮頂罐儲油溫降特點[J].油氣儲運,2003,22(9):47-49.YU D,FANG X Y,LI D F,et al.Temperature drop rule for the large scaled floating roof oil tank[J].Oil & Gas Storage and Transportation,2003,22(9):47-49.

    [26]于達(dá).大型浮頂油罐測溫系統(tǒng)的研發(fā)[J].油氣儲運,2005,24(8):41-43.YU D.Development of the temperature test system for the large scaled floating roof oil tank[J].Oil & Gas Storage and Transportation,2005,24(8):41-43.

    date:2017-06-19.

    ZHAO Jian,soulkissing@163.com

    supported by the National Natural Science Foundation of China (51674086,51704077),the University Nursing Program for Young Scholars with Creative Talents in Heilongjiang Province of China(UNPYSCT-2016125)and the Fostering Foundation for the University Nursing Program for Young Scholars with Creative Talents in Heilongjiang Province of China (SCXHB201601).

    Cooling and gelatinization rule of waxy crude oil during static storage in floating roof tank

    ZHAO Jian1,DONG Hang1,FU Xiaoming2,LEI Qimeng1
    (1Key Laboratory for Enhanced Oil and Gas Recovery of Northeast Petroleum University,Daqing163318,Heilongjiang,China;2Transportation and Storage and Marketing Subsidiary Company of Daqing Oil Field,Daqing163000,Heilongjiang,China)

    Based on the assumption that waxy crude oil is deemed as the porous medium system under the temperature of losing flow point,the additional specific heat capacity method is used to describe the wax crystallization heat,and the momentum source term method is adopted to represent the flow resistance on the liquid oil due to the formation of network structure made by waxy crystal.The finite volume method is used to simulate the gelatinization process of waxy crude oil.The result shows that the gelatinization evolution is mainly conducted by the heat transfer mechanism and boundary condition.Affected by natural convection,the gelled crude oil first generates in the region surrounded by the sidewall and base wall,which is always the most serious gelled region in the tank.The integrated gelled oil layer first appears on the top wall,then it experiences two steps which include the slow growth and fast growth.Moreover,the gelled oil layer is gradually tends to be even when the cooling advances.Later after it appears on the top wall,the gelled oil layer covers the base wall which has the opposite development process from that on the top wall.The gelled oil layer covers the sidewall at the last step.On the sidewall the evolution of gelled oil layer has the trait of advancing from the base wall to the top wall along the sidewall.When natural convection gets stronger,the increasing rate of gelled oil layer turns to be smaller on the top wall,but it has the opposite behavior on the base wall.Based on the temperature profile and evolution of gelled structure,the cooling process can be divided into three states.Natural convection dominates the heat transfer behavior in the first stage.In the second stage,heat conduction gradually replaced the dominant role of natural convection.Finally the heat transfer behavior is conducted by heat conduction and boundary condition in the third stage.In the meantime,the rule of temperature profile and heat loss is also provided.

    crude oil; heat transfer; static storage; gelatinization; numerical simulation

    TE 832.2

    A

    0438—1157(2017)12—4882—10

    10.11949/j.issn.0438-1157.20170788

    2017-06-19收到初稿,2017-09-24收到修改稿。

    聯(lián)系人及第一作者:趙健(1986—),男,博士,講師。

    國家自然科學(xué)基金項目(51674086,51704077);黑龍江省普通本科高等學(xué)校青年創(chuàng)新人才培養(yǎng)計劃項目(UNPYSCT-2016125);東北石油大學(xué)優(yōu)秀科研人才培育基金項目 (SCXHB201601)。

    猜你喜歡
    含蠟溫降罐底
    太古長輸供熱管線溫降統(tǒng)計及分析
    高含蠟原油管道蠟沉積研究綜述
    云南化工(2021年5期)2021-12-21 07:41:14
    汽輪機(jī)快冷新方式在實踐中的成功應(yīng)用
    含蠟原油靜態(tài)磁處理降黏試驗
    化工管理(2020年20期)2020-07-25 15:02:34
    臨邑至濟(jì)南原油管道清管效果研究
    石油研究(2019年5期)2019-09-10 16:17:24
    天然氣與石油(2018年4期)2018-09-08 05:51:06
    集輸管道效率分析及相關(guān)計算探討
    揚州石化儲罐罐底腐蝕及防治
    大型立式儲油罐雙層鋼制罐底結(jié)構(gòu)設(shè)計探討
    化工管理(2017年9期)2017-03-05 12:05:20
    含蠟原油非牛頓流變特性
    精品久久久久久电影网| 一级毛片黄色毛片免费观看视频| 18禁动态无遮挡网站| 看十八女毛片水多多多| 国产 精品1| 免费av中文字幕在线| 国产成人精品久久久久久| 精品国产一区二区三区四区第35| 老司机深夜福利视频在线观看 | 国产麻豆69| 国产精品国产av在线观看| 黑人欧美特级aaaaaa片| 麻豆av在线久日| h视频一区二区三区| 大片免费播放器 马上看| 亚洲欧美成人综合另类久久久| 亚洲成人国产一区在线观看 | √禁漫天堂资源中文www| 国产男女超爽视频在线观看| 在线观看www视频免费| 搡老乐熟女国产| 交换朋友夫妻互换小说| 午夜福利视频在线观看免费| 好男人视频免费观看在线| 久久这里只有精品19| 久久久国产精品麻豆| 在线免费观看不下载黄p国产| 中国国产av一级| a级毛片在线看网站| 在线观看一区二区三区激情| 曰老女人黄片| 久久ye,这里只有精品| 又大又爽又粗| 亚洲精品一二三| 97精品久久久久久久久久精品| 久久影院123| 欧美人与性动交α欧美精品济南到| 美女高潮到喷水免费观看| 青春草国产在线视频| 观看美女的网站| 少妇人妻久久综合中文| 久久精品国产亚洲av高清一级| 国产淫语在线视频| 亚洲四区av| 欧美精品一区二区免费开放| 黄片小视频在线播放| 青春草国产在线视频| 观看美女的网站| 日本一区二区免费在线视频| 亚洲人成77777在线视频| 十八禁高潮呻吟视频| 嫩草影院入口| 国产男人的电影天堂91| 满18在线观看网站| 热99久久久久精品小说推荐| 国产精品熟女久久久久浪| 男女高潮啪啪啪动态图| 精品国产乱码久久久久久小说| 免费高清在线观看日韩| 精品一品国产午夜福利视频| 久热爱精品视频在线9| 中文字幕最新亚洲高清| 又大又爽又粗| 高清视频免费观看一区二区| 国产日韩欧美视频二区| 91老司机精品| 亚洲av国产av综合av卡| 成人免费观看视频高清| 亚洲美女搞黄在线观看| 久久精品国产综合久久久| 夫妻午夜视频| 国产熟女午夜一区二区三区| 亚洲av电影在线进入| 国产乱人偷精品视频| 夜夜骑夜夜射夜夜干| 国产激情久久老熟女| 一级,二级,三级黄色视频| 精品人妻在线不人妻| 精品国产超薄肉色丝袜足j| 黄色毛片三级朝国网站| 最近最新中文字幕免费大全7| 香蕉国产在线看| 秋霞伦理黄片| 亚洲av成人不卡在线观看播放网 | 国产亚洲午夜精品一区二区久久| 男女午夜视频在线观看| 一区在线观看完整版| 欧美最新免费一区二区三区| 国产xxxxx性猛交| 性色av一级| 日本爱情动作片www.在线观看| 日韩不卡一区二区三区视频在线| 国产av精品麻豆| 少妇精品久久久久久久| 下体分泌物呈黄色| 国精品久久久久久国模美| 久久这里只有精品19| 亚洲国产精品国产精品| 欧美久久黑人一区二区| 在线观看免费午夜福利视频| www.av在线官网国产| 国产日韩欧美亚洲二区| 久久久久网色| 国产精品女同一区二区软件| 街头女战士在线观看网站| 亚洲精品美女久久av网站| 成年av动漫网址| www日本在线高清视频| 丰满迷人的少妇在线观看| 久久久久精品久久久久真实原创| 成年美女黄网站色视频大全免费| 国产精品久久久人人做人人爽| 久久久久精品性色| 国产1区2区3区精品| 国产精品三级大全| 在线观看免费午夜福利视频| 亚洲欧美精品自产自拍| 女人久久www免费人成看片| 亚洲第一青青草原| 伊人久久国产一区二区| 亚洲成人一二三区av| 精品国产超薄肉色丝袜足j| av不卡在线播放| 日日爽夜夜爽网站| av网站在线播放免费| 欧美激情高清一区二区三区 | 视频区图区小说| 观看av在线不卡| 国产有黄有色有爽视频| 午夜福利,免费看| 国产午夜精品一二区理论片| 国产亚洲精品第一综合不卡| 18禁动态无遮挡网站| 秋霞在线观看毛片| 亚洲欧美一区二区三区久久| 国产亚洲午夜精品一区二区久久| 菩萨蛮人人尽说江南好唐韦庄| 色精品久久人妻99蜜桃| videosex国产| 99国产精品免费福利视频| 两个人免费观看高清视频| 国产又爽黄色视频| 亚洲人成网站在线观看播放| 欧美激情 高清一区二区三区| 宅男免费午夜| 亚洲精品美女久久久久99蜜臀 | 在线观看免费午夜福利视频| 国产女主播在线喷水免费视频网站| 香蕉丝袜av| 久久久欧美国产精品| 自线自在国产av| 亚洲欧美一区二区三区国产| 最近中文字幕高清免费大全6| av网站在线播放免费| 日韩精品有码人妻一区| 亚洲av日韩精品久久久久久密 | 两个人看的免费小视频| 亚洲av成人精品一二三区| 只有这里有精品99| 狠狠婷婷综合久久久久久88av| 九草在线视频观看| 91国产中文字幕| 18禁国产床啪视频网站| 中文字幕人妻丝袜一区二区 | 色94色欧美一区二区| 一区二区日韩欧美中文字幕| 国产亚洲av片在线观看秒播厂| 国产免费一区二区三区四区乱码| 91老司机精品| 人人澡人人妻人| 黄色怎么调成土黄色| 午夜激情久久久久久久| 国产一区二区三区av在线| 老鸭窝网址在线观看| 国产成人精品福利久久| 在线观看免费视频网站a站| 美女国产高潮福利片在线看| 亚洲,一卡二卡三卡| 精品一区在线观看国产| 色吧在线观看| 日韩精品有码人妻一区| 一级,二级,三级黄色视频| 色网站视频免费| 午夜福利视频在线观看免费| 亚洲成人手机| 成人国语在线视频| 麻豆av在线久日| 午夜福利网站1000一区二区三区| 欧美日本中文国产一区发布| 宅男免费午夜| 亚洲国产看品久久| 亚洲国产精品国产精品| 亚洲精品成人av观看孕妇| 中文字幕色久视频| 这个男人来自地球电影免费观看 | 亚洲av综合色区一区| 亚洲精品国产一区二区精华液| 黄片无遮挡物在线观看| 亚洲精品视频女| av网站免费在线观看视频| 人人妻人人澡人人爽人人夜夜| 欧美乱码精品一区二区三区| 日本一区二区免费在线视频| 婷婷色综合大香蕉| 天美传媒精品一区二区| 亚洲一级一片aⅴ在线观看| 国产一区二区在线观看av| 国产精品秋霞免费鲁丝片| 亚洲av成人精品一二三区| 哪个播放器可以免费观看大片| 成人手机av| 免费人妻精品一区二区三区视频| 18在线观看网站| 精品久久久久久电影网| 精品久久蜜臀av无| 国产亚洲午夜精品一区二区久久| 国产日韩一区二区三区精品不卡| 丰满乱子伦码专区| 菩萨蛮人人尽说江南好唐韦庄| 制服诱惑二区| 亚洲国产欧美一区二区综合| 久久女婷五月综合色啪小说| 久久久国产精品麻豆| 深夜精品福利| 宅男免费午夜| 一区福利在线观看| 天天躁狠狠躁夜夜躁狠狠躁| 国产免费又黄又爽又色| www.自偷自拍.com| 丝瓜视频免费看黄片| 久久精品国产a三级三级三级| 国产亚洲一区二区精品| 欧美成人午夜精品| 纯流量卡能插随身wifi吗| 天天躁夜夜躁狠狠躁躁| 男人舔女人的私密视频| 男男h啪啪无遮挡| 国产高清国产精品国产三级| 黑人欧美特级aaaaaa片| 久久精品国产综合久久久| 久久精品国产综合久久久| 毛片一级片免费看久久久久| av卡一久久| 看免费av毛片| 国产1区2区3区精品| av电影中文网址| 国产一区二区 视频在线| 国产男人的电影天堂91| 黑人猛操日本美女一级片| 免费av中文字幕在线| 少妇 在线观看| 色视频在线一区二区三区| 国产精品嫩草影院av在线观看| 亚洲熟女毛片儿| 97精品久久久久久久久久精品| 亚洲综合精品二区| 亚洲欧美清纯卡通| 日本爱情动作片www.在线观看| 久久久久久久久免费视频了| 久久99热这里只频精品6学生| 久久精品久久精品一区二区三区| 91精品伊人久久大香线蕉| 色播在线永久视频| a 毛片基地| 青春草国产在线视频| 国产无遮挡羞羞视频在线观看| 婷婷色av中文字幕| 啦啦啦啦在线视频资源| 午夜福利影视在线免费观看| av电影中文网址| 一边摸一边抽搐一进一出视频| 一级a爱视频在线免费观看| 亚洲熟女毛片儿| 丝袜在线中文字幕| 亚洲av在线观看美女高潮| 国产日韩欧美亚洲二区| 十八禁高潮呻吟视频| 国产成人系列免费观看| h视频一区二区三区| 免费高清在线观看日韩| 深夜精品福利| 亚洲成av片中文字幕在线观看| 欧美黑人欧美精品刺激| 精品酒店卫生间| 免费高清在线观看日韩| 欧美激情高清一区二区三区 | 丰满乱子伦码专区| 成人影院久久| 精品一区二区三区av网在线观看 | 性少妇av在线| 国产成人91sexporn| 日韩大片免费观看网站| 国产视频首页在线观看| 在线观看免费高清a一片| 天天躁日日躁夜夜躁夜夜| 五月天丁香电影| 999精品在线视频| 乱人伦中国视频| 午夜福利乱码中文字幕| 一本一本久久a久久精品综合妖精| 欧美老熟妇乱子伦牲交| 精品酒店卫生间| 国产精品 国内视频| 一区在线观看完整版| 性少妇av在线| 欧美少妇被猛烈插入视频| 成人影院久久| 9191精品国产免费久久| 女性被躁到高潮视频| 麻豆乱淫一区二区| 女人被躁到高潮嗷嗷叫费观| 九九爱精品视频在线观看| 天天影视国产精品| 91成人精品电影| 国产午夜精品一二区理论片| 黑丝袜美女国产一区| 亚洲精品aⅴ在线观看| 欧美精品一区二区免费开放| 最新在线观看一区二区三区 | av在线app专区| 赤兔流量卡办理| 国产日韩欧美在线精品| 熟女av电影| 99re6热这里在线精品视频| videosex国产| 久久精品国产亚洲av高清一级| 精品免费久久久久久久清纯 | 熟妇人妻不卡中文字幕| 极品少妇高潮喷水抽搐| 国产av一区二区精品久久| 女性生殖器流出的白浆| av网站免费在线观看视频| 亚洲综合色网址| 一本大道久久a久久精品| 久久久国产欧美日韩av| 久久女婷五月综合色啪小说| 精品国产一区二区久久| 精品少妇内射三级| 久久久久久久大尺度免费视频| 人人妻人人爽人人添夜夜欢视频| 亚洲精品aⅴ在线观看| 亚洲精品乱久久久久久| 女人爽到高潮嗷嗷叫在线视频| 晚上一个人看的免费电影| 亚洲国产精品999| 精品免费久久久久久久清纯 | 久久性视频一级片| 十八禁人妻一区二区| 国产免费福利视频在线观看| 免费观看a级毛片全部| 成人三级做爰电影| 嫩草影院入口| 久久久精品区二区三区| 欧美在线黄色| 久久女婷五月综合色啪小说| 欧美在线一区亚洲| 日本vs欧美在线观看视频| 丰满乱子伦码专区| 咕卡用的链子| 久久精品久久精品一区二区三区| 一二三四在线观看免费中文在| 久久精品国产a三级三级三级| 亚洲精品一二三| 天堂中文最新版在线下载| 亚洲色图 男人天堂 中文字幕| √禁漫天堂资源中文www| 91精品伊人久久大香线蕉| 国产精品国产三级国产专区5o| 天天影视国产精品| 波野结衣二区三区在线| 波多野结衣一区麻豆| √禁漫天堂资源中文www| 人妻人人澡人人爽人人| 飞空精品影院首页| 巨乳人妻的诱惑在线观看| 亚洲国产看品久久| 国产精品二区激情视频| 日韩一本色道免费dvd| av片东京热男人的天堂| 男女高潮啪啪啪动态图| 99久久99久久久精品蜜桃| 激情视频va一区二区三区| 日韩 欧美 亚洲 中文字幕| 国产毛片在线视频| 纵有疾风起免费观看全集完整版| 一区在线观看完整版| 男女高潮啪啪啪动态图| 黄色毛片三级朝国网站| 亚洲欧美清纯卡通| 婷婷色综合www| 精品国产国语对白av| 国产亚洲最大av| 欧美国产精品一级二级三级| 成人影院久久| 两个人看的免费小视频| 日本vs欧美在线观看视频| 国产 一区精品| 制服诱惑二区| 免费人妻精品一区二区三区视频| 波野结衣二区三区在线| 久久国产精品大桥未久av| 一区二区日韩欧美中文字幕| 午夜福利乱码中文字幕| 一区二区日韩欧美中文字幕| 国产xxxxx性猛交| 国产色婷婷99| 一级爰片在线观看| 免费黄频网站在线观看国产| 伊人久久国产一区二区| 丁香六月天网| 巨乳人妻的诱惑在线观看| 精品亚洲成国产av| av在线app专区| 亚洲成人免费av在线播放| 人人妻人人爽人人添夜夜欢视频| 国产成人精品无人区| 99久久精品国产亚洲精品| 天美传媒精品一区二区| 国产乱来视频区| 亚洲精品中文字幕在线视频| 搡老岳熟女国产| 人体艺术视频欧美日本| 日本一区二区免费在线视频| 在线观看国产h片| 狂野欧美激情性bbbbbb| 日本91视频免费播放| 老司机影院毛片| 一区在线观看完整版| 人人澡人人妻人| 久热爱精品视频在线9| 亚洲成人免费av在线播放| 国产av精品麻豆| 亚洲美女视频黄频| 亚洲精品av麻豆狂野| 另类精品久久| 亚洲人成网站在线观看播放| 欧美黄色片欧美黄色片| 丝袜人妻中文字幕| 欧美中文综合在线视频| 午夜91福利影院| 老熟女久久久| 亚洲人成网站在线观看播放| 欧美精品人与动牲交sv欧美| 菩萨蛮人人尽说江南好唐韦庄| 好男人视频免费观看在线| 如何舔出高潮| 国产精品二区激情视频| 国产极品天堂在线| 在线观看三级黄色| 久久精品亚洲熟妇少妇任你| 99久久人妻综合| 亚洲,欧美精品.| 伊人亚洲综合成人网| 又大又黄又爽视频免费| 国产成人欧美在线观看 | 成人国语在线视频| 免费av中文字幕在线| 亚洲精品国产av蜜桃| 亚洲欧美精品自产自拍| 午夜久久久在线观看| 久久久久久久大尺度免费视频| 免费观看av网站的网址| 天堂8中文在线网| 国产精品麻豆人妻色哟哟久久| 最新的欧美精品一区二区| 嫩草影院入口| 日韩 欧美 亚洲 中文字幕| 观看av在线不卡| 下体分泌物呈黄色| 日韩,欧美,国产一区二区三区| 爱豆传媒免费全集在线观看| 国产亚洲av高清不卡| www.av在线官网国产| 18禁裸乳无遮挡动漫免费视频| 精品国产国语对白av| 中文字幕亚洲精品专区| 亚洲男人天堂网一区| 纵有疾风起免费观看全集完整版| 十分钟在线观看高清视频www| 91精品伊人久久大香线蕉| 国产欧美日韩一区二区三区在线| 亚洲男人天堂网一区| 亚洲一区二区三区欧美精品| 欧美日韩亚洲国产一区二区在线观看 | 伦理电影免费视频| 9热在线视频观看99| 午夜福利乱码中文字幕| 伊人亚洲综合成人网| 久久99热这里只频精品6学生| 午夜福利乱码中文字幕| 日韩一卡2卡3卡4卡2021年| 黄网站色视频无遮挡免费观看| 免费观看人在逋| svipshipincom国产片| 亚洲欧美日韩另类电影网站| 亚洲精品国产av蜜桃| 丝袜喷水一区| 巨乳人妻的诱惑在线观看| 涩涩av久久男人的天堂| av视频免费观看在线观看| 成人午夜精彩视频在线观看| 一级,二级,三级黄色视频| 好男人视频免费观看在线| 赤兔流量卡办理| 精品国产一区二区久久| 丰满饥渴人妻一区二区三| 日韩熟女老妇一区二区性免费视频| 亚洲欧洲精品一区二区精品久久久 | 天堂俺去俺来也www色官网| 大话2 男鬼变身卡| 国产成人免费无遮挡视频| 亚洲欧美成人精品一区二区| 啦啦啦在线观看免费高清www| 18禁观看日本| 亚洲av欧美aⅴ国产| 一本久久精品| 国产精品.久久久| 亚洲国产精品一区三区| 秋霞伦理黄片| 亚洲av福利一区| 国产在线一区二区三区精| av免费观看日本| 老司机在亚洲福利影院| 欧美日韩综合久久久久久| 午夜免费鲁丝| 国产成人午夜福利电影在线观看| 中文字幕色久视频| 精品人妻一区二区三区麻豆| 如日韩欧美国产精品一区二区三区| 在线 av 中文字幕| 欧美日韩亚洲综合一区二区三区_| 啦啦啦视频在线资源免费观看| 18禁国产床啪视频网站| 在线观看三级黄色| 老熟女久久久| 人人妻人人澡人人看| 如何舔出高潮| 国产精品成人在线| 9热在线视频观看99| 国产成人精品久久久久久| 女人爽到高潮嗷嗷叫在线视频| 啦啦啦 在线观看视频| 伦理电影免费视频| 午夜福利乱码中文字幕| 少妇的丰满在线观看| 十八禁高潮呻吟视频| 国产国语露脸激情在线看| 国产精品蜜桃在线观看| 欧美日韩成人在线一区二区| 2021少妇久久久久久久久久久| 大陆偷拍与自拍| 欧美日韩视频高清一区二区三区二| a级片在线免费高清观看视频| www.精华液| 一级爰片在线观看| 国产精品人妻久久久影院| 国产一区有黄有色的免费视频| 性少妇av在线| 精品人妻一区二区三区麻豆| 日韩制服骚丝袜av| 女人被躁到高潮嗷嗷叫费观| 制服诱惑二区| 久久久精品区二区三区| 国产精品国产三级专区第一集| 亚洲五月色婷婷综合| 久久久亚洲精品成人影院| 水蜜桃什么品种好| 51午夜福利影视在线观看| 亚洲国产最新在线播放| 日韩熟女老妇一区二区性免费视频| 国精品久久久久久国模美| 深夜精品福利| 夫妻午夜视频| 一区二区日韩欧美中文字幕| 成人三级做爰电影| 久久狼人影院| 国产成人啪精品午夜网站| 欧美97在线视频| 高清av免费在线| 国产xxxxx性猛交| 亚洲精品国产一区二区精华液| 亚洲国产看品久久| www.精华液| 性色av一级| 青草久久国产| 老汉色∧v一级毛片| 欧美人与性动交α欧美精品济南到| 久久久久久免费高清国产稀缺| 欧美日韩福利视频一区二区| 熟妇人妻不卡中文字幕| 亚洲一级一片aⅴ在线观看| 国产乱人偷精品视频| 99热网站在线观看| 国产亚洲av片在线观看秒播厂| 一本—道久久a久久精品蜜桃钙片| 亚洲第一av免费看| 免费久久久久久久精品成人欧美视频| 亚洲精品久久成人aⅴ小说| 多毛熟女@视频| 18禁国产床啪视频网站| 精品第一国产精品| 精品人妻熟女毛片av久久网站| 免费观看性生交大片5| 黄色视频不卡| 99九九在线精品视频| 国产av精品麻豆| 婷婷色综合www| 亚洲第一区二区三区不卡| 曰老女人黄片| 又大又爽又粗| 欧美亚洲日本最大视频资源| 国产国语露脸激情在线看| 十八禁高潮呻吟视频| 十分钟在线观看高清视频www| 男女午夜视频在线观看|