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

    多孔介質(zhì)內(nèi)氣泡Ostwald 熟化特性三維孔網(wǎng)數(shù)值模擬*

    2023-09-06 02:15:28張沐安王進(jìn)卿吳睿馮致詹明秀徐旭池作和
    物理學(xué)報(bào) 2023年16期
    關(guān)鍵詞:傳質(zhì)毛細(xì)氣泡

    張沐安 王進(jìn)卿? 吳睿 馮致 詹明秀 徐旭 池作和

    1) (中國計(jì)量大學(xué)計(jì)量測試工程學(xué)院,杭州 310018)

    2) (上海交通大學(xué)機(jī)械與動(dòng)力工程學(xué)院,上海 200240)

    多孔介質(zhì)內(nèi)氣泡的Ostwald 熟化行為廣泛存在于CO2 地質(zhì)封存、多孔材料制備、燃料電池等領(lǐng)域.為探究孔隙尺度下多孔介質(zhì)內(nèi)氣泡的熟化特性,建立了基于濃度耦合計(jì)算的三維孔隙網(wǎng)絡(luò)模型,該模型考慮了氣泡形態(tài)、多孔介質(zhì)結(jié)構(gòu)以及氣液之間的傳質(zhì),通過求解三維孔隙網(wǎng)絡(luò)內(nèi)各孔體的氣相濃度,得到各氣泡的演化過程,并采用四孔隙結(jié)構(gòu)微流體芯片可視化實(shí)驗(yàn)驗(yàn)證了模型的可靠性.為分析多孔介質(zhì)非均質(zhì)性對氣泡熟化過程影響,構(gòu)建了兩種不同孔隙尺寸的三維孔網(wǎng)結(jié)構(gòu),對兩個(gè)區(qū)域內(nèi)的氣泡熟化過程進(jìn)行了數(shù)值模擬.結(jié)果表明: 氣泡初始分布會對熟化過程產(chǎn)生影響,當(dāng)氣泡非均勻分布時(shí),氣泡從小孔隙區(qū)域傳輸至大孔隙區(qū)域的同時(shí),也會各自向自身區(qū)域的大氣泡區(qū)域傳質(zhì);氣泡初始尺寸的差異,會加速熟化進(jìn)程,使熟化時(shí)間明顯短于均勻分布狀態(tài).孔隙數(shù)的選取對平均毛細(xì)力、飽和度等連續(xù)尺度等效參數(shù)具有顯著影響,孔隙數(shù)增加時(shí)毛細(xì)力與飽和度呈現(xiàn)更具規(guī)律性的非線性變化.該模型的建立可以預(yù)測地質(zhì)封存過程中CO2 的演化過程,為CO2長期封存過程中非均質(zhì)性的影響機(jī)制研究提供指導(dǎo).

    1 引言

    氣泡Ostwald 熟化效應(yīng)是一種在氣液兩相流體系中,受氣泡間不同毛細(xì)壓力的影響,小氣泡通過潤濕相分子擴(kuò)散將自身質(zhì)量傳遞給附近大氣泡的現(xiàn)象.毛細(xì)力決定氣泡間濃度分布,濃度梯度驅(qū)動(dòng)的擴(kuò)散效應(yīng)是傳質(zhì)的主導(dǎo)機(jī)制.該過程普遍存在于自然界及能源、化工、生化和環(huán)境保護(hù)等領(lǐng)域,如超疏水表面自清潔技術(shù)[1,2]、油水分離工藝[3]、巖漿脫氣[4,5]、多孔材料制備[6,7]和CO2地質(zhì)封存[8,9]等.最典型的是在CO2地質(zhì)封存中,這種熟化效應(yīng)會引起CO2氣體在局部區(qū)域的重新聚集,增加了地質(zhì)儲存泄漏的風(fēng)險(xiǎn)[10,11].因此,對多孔介質(zhì)內(nèi)氣泡熟化特性的研究至關(guān)重要.

    Ostwald 熟化效應(yīng)可根據(jù)空間狀態(tài)分為自由流體空間與受限空間(多孔介質(zhì))熟化兩種形式.早在20 世紀(jì)眾多國內(nèi)外學(xué)者已經(jīng)對氣泡在溶液中的熟化特性進(jìn)行研究,提出并改進(jìn)了氣泡簇半徑與時(shí)間的穩(wěn)態(tài)生長理論[12-15].在自由流體中,氣泡熟化可以不受空間限制,呈現(xiàn)標(biāo)準(zhǔn)的Ostwald 熟化現(xiàn)象.而在多孔介質(zhì)中,空間結(jié)構(gòu)限制會增加Ostwald熟化過程的復(fù)雜性.目前,國內(nèi)外學(xué)者已對多孔介質(zhì)環(huán)境下的熟化效應(yīng)進(jìn)行了廣泛的實(shí)驗(yàn)研究.Xu 等[16]在實(shí)驗(yàn)過程中成功捕捉到小氣泡與大氣泡的熟化過程,驗(yàn)證了CO2地質(zhì)封存過程中Ostwald 熟化機(jī)制的存在;Xu 等[17]通過2.5D 玻璃微流體芯片實(shí)驗(yàn)發(fā)現(xiàn)氣泡在均質(zhì)多孔介質(zhì)內(nèi)的熟化受限于孔喉的幾何結(jié)構(gòu),最終氣泡尺寸趨于一致,發(fā)生反Ostwald 熟化效應(yīng).Ostwald 熟化是長時(shí)間尺度的過程,在毫米尺度空間域,CO2達(dá)到毛細(xì)平衡需要數(shù)月至數(shù)年的時(shí)間,而對于米級系統(tǒng),熟化則需要一萬年甚至更久[18].可見實(shí)驗(yàn)僅能對短時(shí)間內(nèi)的熟化效應(yīng)進(jìn)行研究,采用數(shù)值模擬方法對長時(shí)間尺度氣泡熟化過程進(jìn)行預(yù)測是有效手段.

    在Ostwald 熟化數(shù)值模型的建立方面,相關(guān)的研究工作主要集中在孔網(wǎng)模型方面.孔網(wǎng)模型的核心思想是將孔隙結(jié)構(gòu)簡化為孔隙和孔喉的組成[19-21],從而對其中內(nèi)部平衡問題和傳質(zhì)計(jì)算進(jìn)行簡化.與直接模擬法(DNS)相比具有計(jì)算范圍廣、計(jì)算效率高的優(yōu)點(diǎn),能夠較好地反映多孔介質(zhì)內(nèi)孔隙結(jié)構(gòu)對擴(kuò)散過程的影響,揭示傳質(zhì)過程孔隙尺度機(jī)理.Xu 等[22]通過雙氣泡孔網(wǎng)模型證明重力會誘導(dǎo)CO2向上遷移影響多孔介質(zhì)內(nèi)Ostwald 熟化,并放大孔隙尺度得到一維連續(xù)介質(zhì)模型以評估CO2地質(zhì)封存安全性;Xu 等[23]基于實(shí)驗(yàn)建立了能準(zhǔn)確預(yù)測氣泡演變的均質(zhì)孔網(wǎng)模型,并指出多孔介質(zhì)內(nèi)孔喉的結(jié)構(gòu)會影響Ostwald 熟化方向;Mehmani 等[24]通過孔網(wǎng)模型模擬混相氣泡的熟化過程,提出一種完全隱式算法求解多孔介質(zhì)內(nèi)被捕獲氣泡的演化過程,并開發(fā)了一種預(yù)測氣泡穩(wěn)定平衡分布狀態(tài)理論,證明氣泡在非均質(zhì)多孔介質(zhì)中的平衡比均勻介質(zhì)中的平衡要慢得多.以上文獻(xiàn)主要在一維和二維條件下開展數(shù)值研究,而自然界和工業(yè)上的問題多為三維多孔介質(zhì).Chalendar 等[25]建立了一種三維錐形孔喉數(shù)值模型,模擬地質(zhì)封存時(shí)CO2氣泡演化,通過算法尋找每個(gè)有效擴(kuò)散路徑并對每條路徑的傳質(zhì)進(jìn)行疊加,實(shí)現(xiàn)了對簡單連續(xù)氣泡群演化時(shí)間尺度的數(shù)值模擬.但該模型隨著熟化過程中氣泡數(shù)量減少,可搜索的有效擴(kuò)散路徑數(shù)量大幅增加,算法復(fù)雜度將呈指數(shù)增長.

    基于質(zhì)量守恒與氣液傳質(zhì)方程,本文構(gòu)建一種通過耦合計(jì)算每個(gè)孔隙中氣相濃度場來獲得氣泡生長速率的孔隙網(wǎng)絡(luò)數(shù)值模型.該模型考慮氣泡形態(tài)、多孔介質(zhì)結(jié)構(gòu)以及氣液之間的傳質(zhì),并通過四孔隙結(jié)構(gòu)微流體芯片可視化實(shí)驗(yàn)驗(yàn)證模型的可靠性.最后以CO2地質(zhì)封存為研究背景,采用該模型對包含兩種不同孔隙尺度多孔介質(zhì)內(nèi)的氣泡熟化過程進(jìn)行三維數(shù)值模擬,分析氣泡初始分布和孔隙數(shù)對熟化的影響特性,并探究孔隙尺度平均毛細(xì)力與飽和度的關(guān)系,為后續(xù)的連續(xù)尺度熟化機(jī)制研究提供指導(dǎo).

    2 數(shù)值計(jì)算模型

    2.1 基本方程

    氣泡熟化過程主要受3 個(gè)方程控制: 拉普拉斯方程、亨利定律和菲克定律:

    式中,pn為氣泡內(nèi)部毛細(xì)壓力,ps為溶液壓力,Ra和Rb分別為氣泡的平面和層面曲率半徑,C為溶質(zhì)的濃度,A為兩氣泡間擴(kuò)散面積,L為兩氣泡間擴(kuò)散距離.本文計(jì)算時(shí)選用的擴(kuò)散面積與擴(kuò)散距離如圖1 所示,其中A1為氣泡1 與喉道切線對應(yīng)的弧面面積,A2為喉道截面積,A3為氣泡2 與喉道切線對應(yīng)的弧面面積;L1為氣泡1 到孔體的距離,L2為喉道長度,L3為氣泡2 到孔體的距離.因此,擴(kuò)散面積A與擴(kuò)散距離L之比為

    圖1 雙氣泡系統(tǒng)傳質(zhì)過程Fig.1.Two-bubble system mass transfer process.

    2.2 計(jì)算模型

    2.2.1 模型框架

    建立的孔網(wǎng)模型假設(shè)熟化過程是準(zhǔn)靜態(tài)的,只考慮氣體擴(kuò)散作用.在整個(gè)孔隙網(wǎng)絡(luò)中,孔體為圓柱體結(jié)構(gòu),喉道為長方體結(jié)構(gòu),當(dāng)氣泡直徑小于孔體深度時(shí),氣泡為球形氣泡,其余情況氣泡為圓柱形氣泡.同時(shí)做出以下假設(shè): 氣泡與液相均處于局部熱力學(xué)平衡狀態(tài),即溶解與生長都是瞬時(shí)的;一個(gè)孔體內(nèi)只能存在一個(gè)氣泡且不可移動(dòng),在整個(gè)溶解或生長過程中,氣泡質(zhì)量中心始終與孔體中心重合;當(dāng)氣泡尺寸達(dá)到孔體尺寸后,氣泡即停止生長,不侵入喉道.

    如前所述,國內(nèi)外學(xué)者已對Ostwald 熟化問題開展了數(shù)值研究,不同于他們建立的數(shù)值模型,本模型通過計(jì)算孔體內(nèi)的氣相濃度分布得到氣泡生長速率,使整個(gè)計(jì)算過程更加簡便,提高了計(jì)算效率.如圖2 所示,該模型會依次判斷各孔體內(nèi)是否存在氣泡以及該孔體周圍充滿液相的孔體個(gè)數(shù),并根據(jù)孔體內(nèi)是否存在氣泡將(3)式中的濃度計(jì)算分為兩種情況: 當(dāng)孔體內(nèi)為氣泡時(shí),根據(jù)拉普拉斯方程和亨利定律求解氣泡濃度;當(dāng)孔體內(nèi)為液相時(shí),基于質(zhì)量守恒方程建立線性方程組求解孔體內(nèi)氣相濃度.以下對兩種情況的濃度計(jì)算進(jìn)行介紹.

    圖2 計(jì)算程序算法流程圖Fig.2.Flow chart for the calculation algorithm.

    2.2.2 孔體內(nèi)為氣泡

    如圖3(a)所示,孔體i內(nèi)存在氣泡i,其相鄰六個(gè)方向的孔體j內(nèi)均存在氣泡j.根據(jù)這些氣泡的初始半徑,代入(1)式和(2)式可以分別得到i和j孔體中氣相濃度Ci與Cj,即

    圖3 數(shù)值計(jì)算模型 (a)孔體i 內(nèi)為氣泡;(b)孔體i 內(nèi)為液相Fig.3.Numerical calculation model: (a) Pore body i with a bubble;(b) pore body i filled with liquid phase.

    2.2.3 孔體內(nèi)為液相

    當(dāng)孔體i內(nèi)充滿液相時(shí),無法依據(jù)拉普拉斯方程與亨利定律求解氣相濃度Ci,需基于質(zhì)量守恒與氣液傳質(zhì)方程建立線性方程組求解.構(gòu)建的孔網(wǎng)模型考慮了孔體連接方式的結(jié)構(gòu)化,即每個(gè)孔體周圍都有6 個(gè)孔隙連接.如圖3(b)所示,孔體i中充滿液相,相鄰6 個(gè)方向的孔體j內(nèi)均存在氣泡j,根據(jù)氣泡j的初始半徑,代入(5)式可分別得到相鄰氣泡的氣相濃度Cj,則充滿液相的孔體i中氣相濃度Ci可基于傳質(zhì)方程(6)求解:

    在整個(gè)三維孔網(wǎng)模型中,孔體內(nèi)氣泡與液相隨機(jī)分布,因此對于孔體i內(nèi)充滿液相的情況,需要依次判斷周圍孔體j是否充滿液相,建立關(guān)于濃度的線性方程組,并采用BiCGSTAB 方法[26]求解出充滿液相的孔體內(nèi)氣相濃度.根據(jù)各孔體的氣相濃度,可以得到氣泡的生長速率.氣泡i通過熟化效應(yīng)在單位時(shí)間變化的質(zhì)量為該氣泡與其相鄰氣泡傳遞質(zhì)量的和,見(7)式.當(dāng)熟化系統(tǒng)中所有孔隙內(nèi)的氣相濃度相等時(shí),熟化結(jié)束,停止計(jì)算.

    2.2.4 時(shí)間步長

    在該程序中,一個(gè)時(shí)間步長為整個(gè)孔隙空間內(nèi)所有氣泡中半徑增大或減小指定長度R'所需要的最短時(shí)間.因此,時(shí)間步長會隨著熟化程序的迭代進(jìn)行自動(dòng)優(yōu)化.在熟化初期,氣泡間尺寸半徑差異大,熟化速率快,時(shí)間步長短;而熟化后期,熟化速率減慢,時(shí)間步長增大.本文計(jì)算模型選取的R'值為0.1 μm.

    (8)式為適用于球形與圓柱形氣泡的體積公式,計(jì)算程序中單次輸出的熟化時(shí)間步長為

    2.3 計(jì)算對象

    在實(shí)際多孔介質(zhì)環(huán)境中,多孔介質(zhì)結(jié)構(gòu)通常呈現(xiàn)非均質(zhì)性,為了探究孔隙結(jié)構(gòu)非均質(zhì)性對熟化過程的影響,選取了由兩種不同孔隙尺度區(qū)域構(gòu)成的三維孔網(wǎng)結(jié)構(gòu)作為計(jì)算對象,對兩個(gè)區(qū)域內(nèi)的氣泡熟化過程進(jìn)行數(shù)值模擬.如圖4 所示,小孔隙區(qū)域、大孔隙區(qū)域孔體半徑分別為20 和50 μm,兩個(gè)區(qū)域深度均為30 μm;所有喉道半徑均為5 μm,深度為5 μm.

    圖4 多孔隙結(jié)構(gòu)熟化模型Fig.4.Multi bubble ripening model.

    在孔隙尺度,氣泡群初始分布會顯著影響熟化過程,氣泡群非均勻分布時(shí)的熟化規(guī)律相較均勻分布有明顯差別[23],因此選取了兩種不同氣泡群初始分布工況(工況1#、2#)進(jìn)行分析(表1).與工況1#相比,工況2#僅改變第3 層氣泡群初始半徑.同時(shí),孔隙數(shù)的選取對平均毛細(xì)力、飽和度等連續(xù)尺度等效參數(shù)的獲取具有顯著影響,為了分析孔隙數(shù)對熟化過程的影響,本文構(gòu)建了孔隙數(shù)分別為325 個(gè)(工況1#)和10800 個(gè)(工況3#)兩種多孔介質(zhì),具體參數(shù)見表1.工況3#氣泡群初始分布與工況1#相同.模擬在常溫常壓條件下進(jìn)行,表2 整理了氣泡熟化模擬時(shí)選用的參數(shù).

    表1 模擬工況參數(shù)Table 1.Simulation parameters of three conditions.

    表2 四孔隙及多孔隙結(jié)構(gòu)氣泡熟化模擬參數(shù)Table 2.Simulation parameters of four bubble and muti bubble ripening system.

    3 實(shí)驗(yàn)方法

    本實(shí)驗(yàn)選用硅基微流體芯片,以去離子水作為液相,CO2作為氣相.芯片分別對硅基底進(jìn)行兩次刻蝕制備出孔隙與喉道圖案,得到具有不同深度流道的2.5D 結(jié)構(gòu),能對芯片內(nèi)的氣泡進(jìn)行夾斷,使每個(gè)孔隙里只有一個(gè)氣泡[17].芯片觀察部分為四孔隙結(jié)構(gòu),具體結(jié)構(gòu)如圖5 所示,其中孔隙深度為30 μm,喉道深度為10 μm.本文采用氣泡熟化可視化實(shí)驗(yàn)裝置對氣泡在四孔隙結(jié)構(gòu)內(nèi)熟化過程進(jìn)行研究,具體的實(shí)驗(yàn)裝置與方法見參考文獻(xiàn)[29].

    圖5 微流控芯片內(nèi)部結(jié)構(gòu)Fig.5.Internal structure of microfluidic chip.

    4 討論部分

    4.1 模型驗(yàn)證

    對四孔隙結(jié)構(gòu)氣泡熟化系統(tǒng)進(jìn)行實(shí)驗(yàn)及模擬研究,結(jié)果如圖6 所示.熟化過程中,雖然氣泡1 的半徑小于氣泡3,但熟化結(jié)束后氣泡3 消失而氣泡1 并未消失.這是由于氣泡1 受到了氣泡2 的質(zhì)量傳輸,而氣泡3 周圍存在半徑更大的氣泡4,使得氣泡3 向氣泡4 傳質(zhì),最終導(dǎo)致上述現(xiàn)象的發(fā)生.可見氣泡的初始分布不同,會導(dǎo)致不同熟化過程的發(fā)生.

    圖6 四孔隙結(jié)構(gòu)氣泡熟化過程 (a)實(shí)驗(yàn)過程;(b)模擬過程Fig.6.Four bubble ripening system: (a) Experiments process;(b) simulations process.

    圖7 為四孔隙結(jié)構(gòu)熟化過程中氣泡半徑變化的模擬和實(shí)驗(yàn)結(jié)果對比,模擬和實(shí)驗(yàn)結(jié)果符合較好.在熟化時(shí)間上,除氣泡3 存在差別,其余氣泡達(dá)到穩(wěn)定的模擬時(shí)間與實(shí)驗(yàn)較為吻合.造成這種差別的原因在于模擬過程中氣泡位置始終位于孔體正中心,而在實(shí)際實(shí)驗(yàn)中,氣泡的位置會與孔體正中心存在一定偏差.隨著熟化的進(jìn)行,模擬和實(shí)驗(yàn)結(jié)果偏差會隨之變大.

    圖7 四孔隙結(jié)構(gòu)氣泡半徑隨時(shí)間變化Fig.7.Evolution of curvature radius in four bubble ripening system.

    4.2 氣泡初始分布對熟化過程影響

    由于工況1#各層孔隙結(jié)構(gòu)和氣泡初始半徑均相同,因此三維空間內(nèi)不同深度方向之間不存在傳質(zhì),故選用單層二維平面數(shù)據(jù)進(jìn)行討論.工況1#中氣泡熟化過程如圖8 所示.觀察圖8 發(fā)現(xiàn):1)由于小孔隙區(qū)域氣泡半徑小,毛細(xì)壓力大,因此在熟化初期小孔隙區(qū)域氣泡半徑顯著減小,氣體從該區(qū)域傳輸至大孔隙區(qū)域;2)熟化過程中,大、小孔隙區(qū)域氣泡的消失順序呈現(xiàn)間隔性,如熟化進(jìn)行至14 h 時(shí),小孔隙區(qū)域殘留氣泡呈現(xiàn)規(guī)則性的間隔分布,4.4×1018h 的大孔隙區(qū)域殘留氣泡分布也呈現(xiàn)相同規(guī)律;3)在熟化末期,小孔區(qū)域殘留氣泡集中在中間區(qū)域,即靠近及遠(yuǎn)離大孔隙區(qū)域的氣泡最先消失,該現(xiàn)象對連續(xù)尺度模型的建立具有指導(dǎo)意義,目前文獻(xiàn)認(rèn)為小孔隙區(qū)域氣泡消失的順序應(yīng)從靠近大孔隙區(qū)域開始逐列向遠(yuǎn)離該區(qū)域的方向推進(jìn).

    圖8 工況1#模擬結(jié)果 (a)氣泡演化過程,其中白色區(qū)域?yàn)镃O2 氣泡,藍(lán)色區(qū)域?yàn)槿ルx子水,灰色區(qū)域?yàn)楣桀w粒;(b) CO2 氣泡曲率半徑變化圖Fig.8.Simulation results of condition 1#: (a) Bubble evolution process,where the white regions are CO2 bubble,blue regions are filled with DI water,and the gray regions are silicon grains;(b) curvature radius variation diagram of each CO2 bubble.

    為了解釋上述現(xiàn)象,采用一維孔隙結(jié)構(gòu)的氣泡熟化過程進(jìn)行說明(圖9).一維孔隙結(jié)構(gòu)由3 個(gè)小孔體和1 個(gè)大孔體構(gòu)成,每個(gè)孔體間通過喉道相連.小孔體內(nèi)的氣泡大小相同,且小于大孔體內(nèi)氣泡.氣泡間傳質(zhì)受毛細(xì)力決定的濃度梯度與擴(kuò)散距離影響.因此,熟化開始時(shí),與大氣泡4 毗鄰的小氣泡3 最先發(fā)生熟化,導(dǎo)致氣泡3 半徑減小,此時(shí)氣泡4、氣泡3 和氣泡2 的半徑大小關(guān)系為R4>R2>R3,故氣泡3 同時(shí)向氣泡4 和氣泡2 傳質(zhì)(見圖9 箭頭方向).氣泡2 獲得傳質(zhì)后半徑增大,使得R2>R1,引起氣泡1 向氣泡2 傳質(zhì).上述傳質(zhì)現(xiàn)象的存在,使得氣泡發(fā)生規(guī)則性的間隔消失.在多孔隙結(jié)構(gòu)(圖8(a))中,這種間隔變化使小孔隙區(qū)域內(nèi)氣泡質(zhì)量往大孔隙區(qū)域傳輸?shù)耐瑫r(shí),還會往中間區(qū)域積累,因此在小孔隙區(qū)域內(nèi)位于中間列的氣泡最慢消失.

    圖9 一維孔隙結(jié)構(gòu)氣泡熟化過程Fig.9.Bubble evolution process of one-dimensional structure.

    工況2#中,由于孔體結(jié)構(gòu)和氣泡初始尺寸的對稱性,第1 層和第5 層的熟化情況相同,第2 層與第4 層的熟化情況相同.工況2#中氣泡熟化過程如圖10 所示.觀察圖10 發(fā)現(xiàn),由于氣泡非均勻分布,熟化過程與工況1#存在以下區(qū)別: 1)工況1#中由于氣泡均勻分布,不存在向不同深度的中間層區(qū)域傳質(zhì)的現(xiàn)象,而工況2#中氣泡從小孔隙區(qū)域傳輸至大孔隙區(qū)域的同時(shí),大、小孔隙區(qū)域也會各自向自身區(qū)域的大氣泡區(qū)域傳質(zhì),因此整個(gè)熟化過程會向中間層區(qū)域逐層推進(jìn);2)如圖10(b),(d)所示,由于氣泡初始分布的差異,小孔隙區(qū)域內(nèi)氣泡熟化進(jìn)行更快,整個(gè)孔隙空間內(nèi)所有氣泡完成熟化的時(shí)間為1.7×1018h,與工況1#中完成熟化的時(shí)間量級上較為相近.進(jìn)一步證明了氣泡熟化是一個(gè)長時(shí)間尺度的過程,CO2地質(zhì)封存中達(dá)到平衡需要數(shù)百萬年時(shí)間[30],時(shí)間尺度的研究可以確定封存過程中CO2的長期儲存能力,評估地質(zhì)封存系統(tǒng)的穩(wěn)定性和長期效果;3)對比工況1# (圖8(b))與工況2# (圖10(b),(d),(f))中氣泡曲率半徑變化發(fā)現(xiàn),由于工況2#中自身熟化現(xiàn)象更顯著,在熟化末期大孔隙區(qū)域內(nèi)大部分氣泡都會充滿孔隙.該模型可以模擬CO2氣泡在地下儲層中的傳輸和分布過程,確定潛在的CO2泄漏路徑和風(fēng)險(xiǎn),為CO2地質(zhì)封存技術(shù)提供重要的指導(dǎo)意義.

    圖10 工況2#模擬結(jié)果 (a)第1 層和第5 層氣泡演化過程;(b)第1 層和第5 層CO2 氣泡曲率半徑變化;(c)第2 層和第4 層氣泡演化過程;(d)第2 層和第4 層CO2 氣泡曲率半徑變化;(e)第3 層氣泡演化過程;(f)第3 層CO2 氣泡曲率半徑變化Fig.10.Simulation results of condition 2#: (a) Bubble evolution process in layers 1 and 5;(b) curvature radius variation of each CO2 bubble in layers 1 and 5;(c) bubble evolution process in layers 2 and 4;(d) curvature radius variation of each CO2 bubble in layers 2 and 4;(e) bubble evolution process in layer 3;(f) curvature radius variation of each CO2 bubble in layer 3.

    4.3 孔隙數(shù)對熟化過程影響

    熟化過程連續(xù)尺度模型參數(shù)的構(gòu)建依賴于孔隙尺度熟化特性的研究,可基于孔隙尺度模擬獲得連續(xù)尺度等效參數(shù),如平均毛細(xì)力、飽和度等參數(shù).孔隙尺度下多孔介質(zhì)孔隙數(shù)的選取對連續(xù)尺度等效參數(shù)的獲取具有顯著影響.為此,本文對工況1#和工況3#熟化過程中CO2氣泡受到的毛細(xì)力與飽和度變化進(jìn)行統(tǒng)計(jì).熟化過程中氣泡受到的平均毛細(xì)力如圖11(a)所示,當(dāng)孔隙數(shù)較少時(shí),個(gè)別氣泡半徑的突變也會對整個(gè)孔隙空間產(chǎn)生較大的影響,使毛細(xì)力曲線走勢振蕩;當(dāng)孔隙數(shù)較多時(shí),毛細(xì)力呈下降趨勢,曲線波動(dòng)較少.熟化過程中飽和度隨時(shí)間的變化如圖11(b)所示,熟化初期小孔隙區(qū)域內(nèi)的氣泡向大孔隙區(qū)域傳質(zhì),導(dǎo)致小孔隙區(qū)域飽和度持續(xù)下降,而在熟化后期大孔隙區(qū)域內(nèi)的氣泡之間也會發(fā)生熟化使大孔隙區(qū)域整體飽和度降低.孔隙數(shù)的增加提供了更多的氣體擴(kuò)散路徑,使大孔隙區(qū)域內(nèi)氣泡自身熟化的情況更復(fù)雜.

    圖11 CO2 氣泡各參數(shù)變化圖 (a)毛細(xì)力;(b) CO2 飽和度Fig.11.Variation diagram of CO2 bubble parameters: (a) Capillary pressure;(b) CO2 saturation.

    在構(gòu)建連續(xù)尺度熟化模型時(shí),毛細(xì)力與飽和度關(guān)系式的選取直接影響熟化模型精度[31].目前國內(nèi)外學(xué)者在建立多孔介質(zhì)內(nèi)CO2熟化連續(xù)尺度模型時(shí),認(rèn)為熟化過程中CO2氣泡毛細(xì)力與飽和度呈線性關(guān)系,如Li 等[18]提出的連續(xù)尺度數(shù)值模型中選用了“Pc=a+bSg”的毛細(xì)力模型.本文分別對兩種工況模擬過程中氣泡毛細(xì)力與飽和度進(jìn)行統(tǒng)計(jì),其變化關(guān)系如圖12 所示,兩種工況下大、小孔隙區(qū)域內(nèi)的毛細(xì)力與飽和度的變化均呈現(xiàn)非線性關(guān)系.對比圖12(a)和圖12(b)發(fā)現(xiàn),當(dāng)孔隙數(shù)增加時(shí),大孔隙區(qū)域和小孔隙區(qū)域內(nèi)毛細(xì)力與飽和度關(guān)系的規(guī)律性顯著提升.這是由于熟化過程中氣泡的變化情況非常復(fù)雜,氣泡所受的毛細(xì)力與飽和度都經(jīng)歷了大幅改變.工況3#中毛細(xì)力與飽和度關(guān)系的變化具有更明顯的規(guī)律性,故對該工況下的數(shù)據(jù)點(diǎn)進(jìn)行非線性擬合,得到小孔隙區(qū)域內(nèi)毛細(xì)力與飽和度的關(guān)系式為大孔隙區(qū)域內(nèi)毛細(xì)力與飽和度的關(guān)系式為Pc=.研究毛細(xì)力和飽和度的關(guān)系,可以提升連續(xù)尺度熟化模型的可靠性和預(yù)測能力,為CO2地質(zhì)封存過程提供指導(dǎo)和決策支持.

    圖12 氣泡毛細(xì)力與CO2 飽和度關(guān)系對比圖 (a)小孔隙區(qū)域;(b)大孔隙區(qū)域Fig.12.Comparison of capillary pressure with CO2 saturation curves: (a) Small pore region;(b) large pore region.

    5 結(jié)論

    本文構(gòu)建了基于濃度耦合計(jì)算的三維孔隙網(wǎng)絡(luò)模型,并采用四孔可視化實(shí)驗(yàn)對模型準(zhǔn)確度進(jìn)行了驗(yàn)證,隨后對三維條件下不同初始分布和不同孔隙數(shù)的氣泡開展數(shù)值模擬研究,結(jié)論如下.

    1)本文建立的模型考慮了氣泡形態(tài)、多孔介質(zhì)結(jié)構(gòu)以及氣液之間的傳質(zhì),通過求解三維孔網(wǎng)內(nèi)各孔體的氣相濃度,可以得到各氣泡的演化過程.該模型數(shù)值模擬結(jié)果與四孔隙結(jié)構(gòu)實(shí)驗(yàn)吻合較好.

    2)數(shù)值模擬結(jié)果表明氣泡初始分布和孔隙數(shù)均會對熟化過程產(chǎn)生影響,當(dāng)氣泡非均勻分布時(shí),由于大、小孔隙區(qū)域會各自向自身區(qū)域的大氣泡傳質(zhì),熟化過程中殘留氣泡不再呈現(xiàn)規(guī)律性的間隔分布,并且熟化進(jìn)程明顯加速;孔隙數(shù)增加時(shí),大、小孔隙區(qū)域的毛細(xì)力與飽和度均呈現(xiàn)更具規(guī)律性的非線性變化.

    3)本文得到大、小孔隙區(qū)域毛細(xì)力與飽和度的關(guān)系式分別為:,與現(xiàn)有文獻(xiàn)的假設(shè)有所不同,該結(jié)果對連續(xù)尺度熟化模型的構(gòu)建具有重要指導(dǎo)意義.但該關(guān)系式并未考慮不同氣泡初始分布、孔隙數(shù)、孔喉比等參數(shù)的影響,后續(xù)將在這些方面開展孔隙尺度模擬.

    猜你喜歡
    傳質(zhì)毛細(xì)氣泡
    “毛細(xì)”政務(wù)號關(guān)停背后
    廉政瞭望(2024年5期)2024-05-26 13:21:07
    檸檬氣泡水
    欣漾(2024年2期)2024-04-27 15:19:49
    SIAU詩杭便攜式氣泡水杯
    新潮電子(2021年7期)2021-08-14 15:53:12
    浮法玻璃氣泡的預(yù)防和控制對策
    多孔建筑材料毛細(xì)吸水系數(shù)連續(xù)測量的方法
    能源工程(2020年6期)2021-01-26 00:55:16
    冰凍氣泡
    出現(xiàn)憋喘 可能是毛細(xì)支氣管炎!
    媽媽寶寶(2017年3期)2017-02-21 01:22:16
    氨基酸鹽吸收二氧化碳過程的傳質(zhì)特性
    高滲鹽水霧化吸入治療毛細(xì)支氣管炎的療效觀察
    PTFE膜吸收CO2的工藝參數(shù)對傳質(zhì)性能的影響
    丰满的人妻完整版| 成年人黄色毛片网站| 日韩 欧美 亚洲 中文字幕| 亚洲午夜精品一区,二区,三区| 国内久久婷婷六月综合欲色啪| 国产精品亚洲av一区麻豆| 久久久久国产一级毛片高清牌| 两人在一起打扑克的视频| 大型av网站在线播放| 色播亚洲综合网| 天堂√8在线中文| 免费电影在线观看免费观看| 丰满人妻一区二区三区视频av | 亚洲熟妇熟女久久| 午夜免费观看网址| 长腿黑丝高跟| 99久久精品热视频| 少妇人妻一区二区三区视频| 日韩欧美国产一区二区入口| 村上凉子中文字幕在线| 全区人妻精品视频| 三级国产精品欧美在线观看 | 国产高清视频在线播放一区| 好男人电影高清在线观看| 天天躁狠狠躁夜夜躁狠狠躁| 草草在线视频免费看| 美女午夜性视频免费| 亚洲av美国av| 99在线人妻在线中文字幕| 国模一区二区三区四区视频 | 国产精品亚洲美女久久久| 亚洲全国av大片| 国产av不卡久久| 桃红色精品国产亚洲av| 国内精品久久久久精免费| 亚洲人成网站在线播放欧美日韩| 亚洲,欧美精品.| 99riav亚洲国产免费| 亚洲国产精品久久男人天堂| 国产精品1区2区在线观看.| 欧美在线黄色| 欧美日韩乱码在线| 欧美日韩亚洲国产一区二区在线观看| 成年版毛片免费区| 一进一出好大好爽视频| 国产亚洲av嫩草精品影院| 午夜成年电影在线免费观看| 亚洲人成网站在线播放欧美日韩| 国产精品久久久久久精品电影| 欧美在线一区亚洲| 丰满人妻熟妇乱又伦精品不卡| 亚洲一区中文字幕在线| 日韩欧美一区二区三区在线观看| 日韩大码丰满熟妇| 一级毛片女人18水好多| 悠悠久久av| www.999成人在线观看| 麻豆av在线久日| 日本在线视频免费播放| 亚洲电影在线观看av| 亚洲成人国产一区在线观看| 亚洲欧美日韩高清专用| 男男h啪啪无遮挡| 999久久久国产精品视频| 久久精品影院6| 精品欧美国产一区二区三| 男女那种视频在线观看| 亚洲精品中文字幕在线视频| 色综合站精品国产| 老司机午夜十八禁免费视频| 91九色精品人成在线观看| 午夜福利在线在线| 露出奶头的视频| 国产成人一区二区三区免费视频网站| 久久精品亚洲精品国产色婷小说| 亚洲精品久久成人aⅴ小说| 男女视频在线观看网站免费 | 久久久久精品国产欧美久久久| 天天一区二区日本电影三级| 国产免费av片在线观看野外av| 国产精品 国内视频| 法律面前人人平等表现在哪些方面| 亚洲午夜精品一区,二区,三区| av欧美777| 免费人成视频x8x8入口观看| 欧美极品一区二区三区四区| 在线十欧美十亚洲十日本专区| 欧美av亚洲av综合av国产av| 男女之事视频高清在线观看| 色综合欧美亚洲国产小说| 欧美 亚洲 国产 日韩一| 天天躁狠狠躁夜夜躁狠狠躁| 1024视频免费在线观看| 国产人伦9x9x在线观看| 成人av一区二区三区在线看| 久久久精品欧美日韩精品| 精品无人区乱码1区二区| 看黄色毛片网站| 久久久久亚洲av毛片大全| or卡值多少钱| 色av中文字幕| 一区二区三区激情视频| 色综合站精品国产| 亚洲欧美精品综合久久99| 日本熟妇午夜| 久久精品国产清高在天天线| 亚洲人成电影免费在线| 成人永久免费在线观看视频| 国产精品久久久av美女十八| 又大又爽又粗| 天堂av国产一区二区熟女人妻 | 日韩欧美免费精品| 亚洲国产精品久久男人天堂| 老司机靠b影院| 88av欧美| e午夜精品久久久久久久| 搞女人的毛片| 国产精品,欧美在线| 久久性视频一级片| 9191精品国产免费久久| 色噜噜av男人的天堂激情| 老鸭窝网址在线观看| 国产精品永久免费网站| 国产精品综合久久久久久久免费| 国产亚洲精品久久久久5区| 久久久久国产一级毛片高清牌| 激情在线观看视频在线高清| 久久久久久久久久黄片| 国内精品久久久久久久电影| 婷婷亚洲欧美| 亚洲av成人av| 99久久精品国产亚洲精品| 亚洲专区字幕在线| 久久精品国产清高在天天线| 中文字幕精品亚洲无线码一区| 9191精品国产免费久久| 日韩欧美一区二区三区在线观看| 可以在线观看的亚洲视频| 色综合婷婷激情| 日韩成人在线观看一区二区三区| 一本久久中文字幕| 久久精品国产综合久久久| 亚洲国产精品999在线| 久久久久久久精品吃奶| 免费观看精品视频网站| 国产成人系列免费观看| 欧美成狂野欧美在线观看| 丁香六月欧美| 男人舔女人下体高潮全视频| 久久香蕉激情| 婷婷亚洲欧美| 国产一区二区激情短视频| 亚洲欧美激情综合另类| 精品少妇一区二区三区视频日本电影| 99在线人妻在线中文字幕| 国内精品久久久久久久电影| 国产精品免费视频内射| 国产一区二区三区在线臀色熟女| 丰满人妻熟妇乱又伦精品不卡| 欧美日韩黄片免| 九九热线精品视视频播放| 首页视频小说图片口味搜索| 亚洲 欧美 日韩 在线 免费| 叶爱在线成人免费视频播放| 亚洲天堂国产精品一区在线| av天堂在线播放| 久久久久久久精品吃奶| 日日摸夜夜添夜夜添小说| 男女午夜视频在线观看| 桃红色精品国产亚洲av| 国产伦在线观看视频一区| 亚洲美女视频黄频| 午夜激情福利司机影院| 国语自产精品视频在线第100页| 日本精品一区二区三区蜜桃| 老司机午夜福利在线观看视频| 国产精品1区2区在线观看.| 免费高清视频大片| 成年女人毛片免费观看观看9| 一二三四在线观看免费中文在| 免费高清视频大片| 国产精品久久电影中文字幕| 好男人电影高清在线观看| 精品免费久久久久久久清纯| 看免费av毛片| 国内揄拍国产精品人妻在线| 午夜老司机福利片| 又紧又爽又黄一区二区| or卡值多少钱| 91老司机精品| 亚洲国产看品久久| bbb黄色大片| 99国产精品一区二区三区| 丰满的人妻完整版| 首页视频小说图片口味搜索| 悠悠久久av| 搡老妇女老女人老熟妇| 成人一区二区视频在线观看| 欧美黑人精品巨大| 亚洲欧美激情综合另类| 999久久久精品免费观看国产| 久久久精品欧美日韩精品| 欧洲精品卡2卡3卡4卡5卡区| 亚洲精品在线观看二区| 香蕉国产在线看| 欧美日本视频| 长腿黑丝高跟| 久久久久精品国产欧美久久久| 亚洲专区中文字幕在线| 中国美女看黄片| 亚洲av中文字字幕乱码综合| 国产精品av视频在线免费观看| 国产真实乱freesex| 欧美 亚洲 国产 日韩一| 日韩欧美三级三区| 亚洲欧美精品综合久久99| 国产一级毛片七仙女欲春2| 免费看美女性在线毛片视频| 国内精品久久久久精免费| 美女大奶头视频| 给我免费播放毛片高清在线观看| 国产伦在线观看视频一区| 美女黄网站色视频| 国产v大片淫在线免费观看| 亚洲精品在线美女| 国产精品 国内视频| 国产精品野战在线观看| 国产精品一及| 人妻久久中文字幕网| 99久久久亚洲精品蜜臀av| 久久精品国产亚洲av高清一级| 久久久久久久久免费视频了| avwww免费| av福利片在线| 亚洲成人免费电影在线观看| 欧美在线一区亚洲| 一本综合久久免费| 久久香蕉国产精品| 久久热在线av| 国产主播在线观看一区二区| 成人国语在线视频| 此物有八面人人有两片| 午夜精品在线福利| av福利片在线观看| 亚洲欧美一区二区三区黑人| 中文字幕av在线有码专区| 国产精品av视频在线免费观看| 久久久精品国产亚洲av高清涩受| 婷婷亚洲欧美| 久久99热这里只有精品18| 亚洲精品国产精品久久久不卡| 51午夜福利影视在线观看| 两人在一起打扑克的视频| 成年女人毛片免费观看观看9| 国产成人一区二区三区免费视频网站| 99久久国产精品久久久| 在线看三级毛片| 久久久久久亚洲精品国产蜜桃av| 国产精品美女特级片免费视频播放器 | 亚洲成a人片在线一区二区| 午夜视频精品福利| 成人特级黄色片久久久久久久| 麻豆成人av在线观看| 日韩精品中文字幕看吧| 国产伦在线观看视频一区| www.精华液| 黄色丝袜av网址大全| 国内精品久久久久精免费| 91九色精品人成在线观看| 色老头精品视频在线观看| 国产亚洲精品久久久久久毛片| 午夜福利高清视频| 五月伊人婷婷丁香| 十八禁网站免费在线| 欧美日韩一级在线毛片| 亚洲av片天天在线观看| 午夜福利18| a级毛片在线看网站| 91大片在线观看| 777久久人妻少妇嫩草av网站| 天天躁狠狠躁夜夜躁狠狠躁| 18禁美女被吸乳视频| 亚洲第一欧美日韩一区二区三区| 一本综合久久免费| 久久精品夜夜夜夜夜久久蜜豆 | 又大又爽又粗| 日韩欧美国产在线观看| 欧美成人午夜精品| 久久中文看片网| 十八禁人妻一区二区| 欧美一级毛片孕妇| 久久香蕉精品热| 国产av不卡久久| av福利片在线| 丁香六月欧美| 国内精品久久久久精免费| 一二三四在线观看免费中文在| 三级国产精品欧美在线观看 | 亚洲欧美日韩高清在线视频| 久久精品亚洲精品国产色婷小说| 午夜福利18| av视频在线观看入口| 亚洲成av人片免费观看| 亚洲成a人片在线一区二区| 国产亚洲精品综合一区在线观看 | 日本五十路高清| 国产一区二区激情短视频| 99热6这里只有精品| 精品欧美一区二区三区在线| av超薄肉色丝袜交足视频| 久久性视频一级片| 亚洲人成网站高清观看| xxx96com| 亚洲精品在线美女| 亚洲熟妇中文字幕五十中出| 午夜免费激情av| 亚洲av中文字字幕乱码综合| 久久久久久人人人人人| 久久香蕉国产精品| av福利片在线| 99热这里只有精品一区 | 国产精品野战在线观看| 悠悠久久av| 欧美日韩乱码在线| 天堂影院成人在线观看| 我的老师免费观看完整版| 美女黄网站色视频| 最近最新免费中文字幕在线| 18禁黄网站禁片午夜丰满| 欧洲精品卡2卡3卡4卡5卡区| 欧美高清成人免费视频www| 日韩欧美 国产精品| 床上黄色一级片| 日韩高清综合在线| www国产在线视频色| 搞女人的毛片| 身体一侧抽搐| 久久久久久国产a免费观看| 少妇裸体淫交视频免费看高清 | 亚洲最大成人中文| 欧美成人午夜精品| 亚洲精品在线观看二区| 999精品在线视频| 欧美绝顶高潮抽搐喷水| 成人av一区二区三区在线看| 少妇粗大呻吟视频| 91成年电影在线观看| 亚洲熟妇熟女久久| av视频在线观看入口| 久久久久久亚洲精品国产蜜桃av| 国产伦一二天堂av在线观看| 国产蜜桃级精品一区二区三区| 性色av乱码一区二区三区2| 男女之事视频高清在线观看| 久久国产精品人妻蜜桃| 此物有八面人人有两片| 久久香蕉国产精品| 午夜福利欧美成人| 51午夜福利影视在线观看| 欧美丝袜亚洲另类 | www.999成人在线观看| 日韩大码丰满熟妇| 国产高清videossex| 欧美三级亚洲精品| 国产亚洲精品第一综合不卡| www日本在线高清视频| 亚洲激情在线av| 香蕉国产在线看| 国产精品99久久99久久久不卡| 久久久久久大精品| 午夜成年电影在线免费观看| 日本黄大片高清| 亚洲狠狠婷婷综合久久图片| 免费看日本二区| 欧美不卡视频在线免费观看 | 国内久久婷婷六月综合欲色啪| 9191精品国产免费久久| 人人妻人人澡欧美一区二区| 欧美性猛交╳xxx乱大交人| 精品无人区乱码1区二区| 女人爽到高潮嗷嗷叫在线视频| av中文乱码字幕在线| 99国产精品一区二区蜜桃av| www.熟女人妻精品国产| 日韩免费av在线播放| 18禁美女被吸乳视频| 观看免费一级毛片| 国产亚洲精品第一综合不卡| 国模一区二区三区四区视频 | 在线观看午夜福利视频| 无人区码免费观看不卡| 最新美女视频免费是黄的| 一级黄色大片毛片| 91字幕亚洲| 久久久精品欧美日韩精品| 丰满人妻熟妇乱又伦精品不卡| av中文乱码字幕在线| 久久久久九九精品影院| 好看av亚洲va欧美ⅴa在| 国产激情偷乱视频一区二区| 精品国产超薄肉色丝袜足j| 天天添夜夜摸| 精品一区二区三区av网在线观看| 在线视频色国产色| 久久香蕉激情| 国产视频一区二区在线看| 日本免费a在线| 99精品在免费线老司机午夜| 国内揄拍国产精品人妻在线| 精品无人区乱码1区二区| 亚洲专区国产一区二区| 黄色视频不卡| 99re在线观看精品视频| 免费看美女性在线毛片视频| www国产在线视频色| 可以在线观看毛片的网站| 国产午夜福利久久久久久| 免费高清视频大片| 成人国语在线视频| 免费无遮挡裸体视频| www.999成人在线观看| 日韩高清综合在线| 91老司机精品| 精品久久久久久,| 中国美女看黄片| 在线观看舔阴道视频| 欧美黄色片欧美黄色片| 成年人黄色毛片网站| 久久伊人香网站| 国产精品乱码一区二三区的特点| 亚洲 国产 在线| 国产一区二区激情短视频| 国产精品亚洲av一区麻豆| 亚洲av成人一区二区三| 在线观看一区二区三区| 欧美人与性动交α欧美精品济南到| 白带黄色成豆腐渣| 精品久久久久久久毛片微露脸| 日本一本二区三区精品| 午夜a级毛片| 国产日本99.免费观看| 国产精品综合久久久久久久免费| 91国产中文字幕| av福利片在线观看| 国产精品自产拍在线观看55亚洲| 日韩欧美在线二视频| 国产爱豆传媒在线观看 | 老鸭窝网址在线观看| 九色国产91popny在线| 18禁观看日本| 中文字幕熟女人妻在线| aaaaa片日本免费| 亚洲色图av天堂| 久久香蕉国产精品| 两人在一起打扑克的视频| 夜夜看夜夜爽夜夜摸| 亚洲性夜色夜夜综合| 久久午夜综合久久蜜桃| 亚洲国产精品999在线| 亚洲成人久久爱视频| 久久久久久免费高清国产稀缺| 99国产精品一区二区三区| 国产亚洲精品一区二区www| 久久天堂一区二区三区四区| 欧美中文日本在线观看视频| 日本黄大片高清| 一边摸一边抽搐一进一小说| 超碰成人久久| 国产午夜精品久久久久久| 制服人妻中文乱码| 久久亚洲精品不卡| 国产v大片淫在线免费观看| 欧美又色又爽又黄视频| 成人av在线播放网站| 日韩欧美 国产精品| tocl精华| 在线观看www视频免费| 一级片免费观看大全| 日本熟妇午夜| 日日干狠狠操夜夜爽| 欧美日韩瑟瑟在线播放| 欧美成人性av电影在线观看| 女人爽到高潮嗷嗷叫在线视频| 国产精品爽爽va在线观看网站| 欧美最黄视频在线播放免费| 国产高清视频在线观看网站| 两个人看的免费小视频| 久久婷婷人人爽人人干人人爱| 日韩欧美在线二视频| 一区福利在线观看| 欧美极品一区二区三区四区| 啦啦啦免费观看视频1| 欧美中文日本在线观看视频| 亚洲美女黄片视频| 后天国语完整版免费观看| 在线a可以看的网站| 亚洲精品美女久久av网站| 久久伊人香网站| 神马国产精品三级电影在线观看 | 成人av在线播放网站| 久久久久久久久中文| 丰满人妻一区二区三区视频av | 两人在一起打扑克的视频| 99在线人妻在线中文字幕| 真人一进一出gif抽搐免费| 久久久久国产精品人妻aⅴ院| 久久婷婷人人爽人人干人人爱| 亚洲第一欧美日韩一区二区三区| 精品久久久久久,| 国产69精品久久久久777片 | 丝袜美腿诱惑在线| 日日爽夜夜爽网站| 草草在线视频免费看| 久久这里只有精品19| 日本成人三级电影网站| 岛国在线免费视频观看| 天堂av国产一区二区熟女人妻 | 丝袜美腿诱惑在线| 免费一级毛片在线播放高清视频| 免费在线观看完整版高清| 日韩欧美在线二视频| 男人舔女人下体高潮全视频| 精品午夜福利视频在线观看一区| 2021天堂中文幕一二区在线观| aaaaa片日本免费| av天堂在线播放| 欧美色欧美亚洲另类二区| 国产高清视频在线播放一区| 欧美色视频一区免费| 色哟哟哟哟哟哟| 午夜精品久久久久久毛片777| 少妇粗大呻吟视频| 在线观看免费午夜福利视频| 12—13女人毛片做爰片一| 香蕉国产在线看| 日本黄大片高清| av福利片在线| 91字幕亚洲| 两人在一起打扑克的视频| 亚洲自拍偷在线| 久久精品国产99精品国产亚洲性色| 久久人人精品亚洲av| e午夜精品久久久久久久| 村上凉子中文字幕在线| 久久精品国产99精品国产亚洲性色| 亚洲成人中文字幕在线播放| 可以免费在线观看a视频的电影网站| 欧美在线黄色| 色尼玛亚洲综合影院| 国产91精品成人一区二区三区| 国产又黄又爽又无遮挡在线| 亚洲熟妇中文字幕五十中出| 成人一区二区视频在线观看| 亚洲国产欧美人成| 麻豆一二三区av精品| 精品久久久久久久毛片微露脸| 黄色毛片三级朝国网站| 床上黄色一级片| 精品第一国产精品| 国产精品一及| 日韩欧美在线二视频| 免费电影在线观看免费观看| 观看免费一级毛片| 国产高清视频在线观看网站| 精品国内亚洲2022精品成人| 亚洲激情在线av| 黑人欧美特级aaaaaa片| 99热只有精品国产| 精品一区二区三区av网在线观看| 丰满人妻一区二区三区视频av | 久热爱精品视频在线9| 男人的好看免费观看在线视频 | 欧美成人一区二区免费高清观看 | 变态另类成人亚洲欧美熟女| 中国美女看黄片| 精品国产超薄肉色丝袜足j| 一本一本综合久久| 精品电影一区二区在线| 日本撒尿小便嘘嘘汇集6| 日韩成人在线观看一区二区三区| 精品日产1卡2卡| 男女下面进入的视频免费午夜| 午夜久久久久精精品| 美女高潮喷水抽搐中文字幕| 欧美一级a爱片免费观看看 | 两人在一起打扑克的视频| 人人妻,人人澡人人爽秒播| 亚洲免费av在线视频| 国产野战对白在线观看| 狂野欧美激情性xxxx| 亚洲国产看品久久| 伦理电影免费视频| 亚洲成人国产一区在线观看| 两个人看的免费小视频| 国产av不卡久久| 啦啦啦观看免费观看视频高清| 国产精品九九99| 黑人操中国人逼视频| 久久久精品欧美日韩精品| 国产精品亚洲美女久久久| 一个人观看的视频www高清免费观看 | 精品免费久久久久久久清纯| 亚洲va日本ⅴa欧美va伊人久久| 搡老妇女老女人老熟妇| 日韩精品中文字幕看吧| 啪啪无遮挡十八禁网站| 久久香蕉精品热| 麻豆成人午夜福利视频| 日本免费一区二区三区高清不卡| 国产成人av教育| 天天躁狠狠躁夜夜躁狠狠躁| 国产激情久久老熟女| 色综合亚洲欧美另类图片| 国产探花在线观看一区二区| or卡值多少钱|