• <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ì)性能的影響
    99热全是精品| 国产精品乱码一区二三区的特点| 2022亚洲国产成人精品| 国产精品日韩av在线免费观看| 一本久久中文字幕| 国产综合懂色| 亚洲人成网站在线播放欧美日韩| 丰满的人妻完整版| 成人无遮挡网站| 欧美成人精品欧美一级黄| 国产精品人妻久久久久久| 国产av在哪里看| 久久久久久久久大av| 亚洲七黄色美女视频| 欧美色欧美亚洲另类二区| 国产伦精品一区二区三区四那| 午夜福利在线观看免费完整高清在 | 日韩中字成人| 久久6这里有精品| 日本黄色片子视频| 久久久久免费精品人妻一区二区| 久久99蜜桃精品久久| 91久久精品国产一区二区三区| 1024手机看黄色片| 美女被艹到高潮喷水动态| 不卡一级毛片| 九色成人免费人妻av| 毛片女人毛片| 亚洲三级黄色毛片| 两个人的视频大全免费| 哪里可以看免费的av片| 久久精品夜色国产| 看片在线看免费视频| 国产精品一区二区性色av| 99热只有精品国产| 欧美色视频一区免费| 久久久a久久爽久久v久久| 国产高清不卡午夜福利| 久久久久久久久久黄片| av女优亚洲男人天堂| 亚洲人成网站高清观看| 男的添女的下面高潮视频| 国产高清视频在线观看网站| 老司机影院成人| 亚洲av成人精品一区久久| 亚洲国产高清在线一区二区三| 亚洲综合色惰| 男人和女人高潮做爰伦理| 日本一本二区三区精品| 男的添女的下面高潮视频| 国产成人福利小说| 欧美色欧美亚洲另类二区| 亚洲成av人片在线播放无| 男人舔奶头视频| 色综合站精品国产| 成人鲁丝片一二三区免费| 精品久久久久久久人妻蜜臀av| 22中文网久久字幕| 舔av片在线| 日韩一区二区三区影片| 国产成人freesex在线| 91av网一区二区| 国产黄片视频在线免费观看| 亚洲人成网站在线播| 黄色一级大片看看| 久久久久久久久大av| 一本久久中文字幕| 麻豆一二三区av精品| 亚洲激情五月婷婷啪啪| 综合色丁香网| 别揉我奶头 嗯啊视频| 国语自产精品视频在线第100页| 一区福利在线观看| 日本一二三区视频观看| a级一级毛片免费在线观看| 男的添女的下面高潮视频| 男女边吃奶边做爰视频| 麻豆av噜噜一区二区三区| 国产午夜精品论理片| 亚洲欧美日韩无卡精品| 色播亚洲综合网| 亚洲美女搞黄在线观看| 婷婷精品国产亚洲av| 精品一区二区免费观看| 久久久色成人| 久久亚洲国产成人精品v| 国产精品无大码| а√天堂www在线а√下载| 91精品国产九色| 国产精品一区二区三区四区免费观看| 菩萨蛮人人尽说江南好唐韦庄 | 国产精品久久久久久精品电影小说 | 中文字幕久久专区| 日本五十路高清| 午夜亚洲福利在线播放| 国产av一区在线观看免费| 国产亚洲精品久久久久久毛片| 精品午夜福利在线看| 在线免费十八禁| 国产在线男女| 国产一区二区在线观看日韩| 国产亚洲精品久久久com| 国产在线精品亚洲第一网站| 一进一出抽搐gif免费好疼| 男女边吃奶边做爰视频| 91精品国产九色| 国语自产精品视频在线第100页| 高清日韩中文字幕在线| 日韩 亚洲 欧美在线| 在线播放无遮挡| av免费观看日本| 男女视频在线观看网站免费| 啦啦啦韩国在线观看视频| 亚洲精品乱码久久久久久按摩| 舔av片在线| 日韩一区二区三区影片| 免费黄网站久久成人精品| 白带黄色成豆腐渣| 久久久午夜欧美精品| 哪里可以看免费的av片| 欧美成人精品欧美一级黄| 夫妻性生交免费视频一级片| 色综合站精品国产| 亚洲美女搞黄在线观看| 亚洲电影在线观看av| 又粗又硬又长又爽又黄的视频 | 全区人妻精品视频| 欧美一区二区国产精品久久精品| 一卡2卡三卡四卡精品乱码亚洲| 欧美区成人在线视频| 亚洲av成人av| 一级av片app| 亚洲av二区三区四区| 亚洲av免费在线观看| 国产黄色视频一区二区在线观看 | 亚洲自拍偷在线| 欧美在线一区亚洲| 99精品在免费线老司机午夜| 国产不卡一卡二| 日韩欧美精品免费久久| kizo精华| 热99在线观看视频| 能在线免费观看的黄片| 国产毛片a区久久久久| 亚洲av男天堂| 国产麻豆成人av免费视频| 久久久久久久久久成人| 爱豆传媒免费全集在线观看| 干丝袜人妻中文字幕| 国产精品乱码一区二三区的特点| 天堂√8在线中文| 国产日本99.免费观看| 黄色配什么色好看| 亚洲内射少妇av| 性欧美人与动物交配| 亚洲色图av天堂| 亚洲精品色激情综合| 中文字幕精品亚洲无线码一区| 欧美人与善性xxx| 久久精品国产亚洲av香蕉五月| 久久国内精品自在自线图片| 欧美+日韩+精品| 日日撸夜夜添| 麻豆av噜噜一区二区三区| 在线观看66精品国产| 99在线人妻在线中文字幕| 久久精品久久久久久久性| 老熟妇乱子伦视频在线观看| 久久久久久大精品| 日韩欧美国产在线观看| 欧美又色又爽又黄视频| 97热精品久久久久久| 国产精品电影一区二区三区| 久久亚洲精品不卡| 少妇人妻精品综合一区二区 | 久久综合国产亚洲精品| 日韩视频在线欧美| 热99在线观看视频| 欧美色欧美亚洲另类二区| 波野结衣二区三区在线| 国产一区二区三区在线臀色熟女| 中文字幕人妻熟人妻熟丝袜美| 欧美一区二区精品小视频在线| 校园人妻丝袜中文字幕| 国产成人精品一,二区 | 最近最新中文字幕大全电影3| 日本免费一区二区三区高清不卡| 午夜福利在线观看吧| 成人综合一区亚洲| 日韩欧美精品v在线| 一本久久精品| 欧美高清成人免费视频www| 少妇丰满av| 日韩亚洲欧美综合| 国产白丝娇喘喷水9色精品| 国产亚洲91精品色在线| 日本三级黄在线观看| 国内少妇人妻偷人精品xxx网站| 校园人妻丝袜中文字幕| 麻豆一二三区av精品| 男女做爰动态图高潮gif福利片| 精品一区二区三区视频在线| 精品不卡国产一区二区三区| 我要看日韩黄色一级片| 精品久久久久久久久久免费视频| 午夜a级毛片| 久99久视频精品免费| 久久鲁丝午夜福利片| 国产精品一及| 欧美3d第一页| 天天躁夜夜躁狠狠久久av| 国产成人a∨麻豆精品| 赤兔流量卡办理| 久久99热这里只有精品18| 51国产日韩欧美| 白带黄色成豆腐渣| 欧美日韩乱码在线| a级一级毛片免费在线观看| 男女做爰动态图高潮gif福利片| 51国产日韩欧美| 欧美最新免费一区二区三区| 午夜亚洲福利在线播放| 在线国产一区二区在线| 亚洲av成人av| 最近最新中文字幕大全电影3| 久久99蜜桃精品久久| 国内精品宾馆在线| 黄色一级大片看看| 一级毛片久久久久久久久女| 看片在线看免费视频| 日韩欧美一区二区三区在线观看| 久久草成人影院| 人妻夜夜爽99麻豆av| 女人十人毛片免费观看3o分钟| 国产精品野战在线观看| 日韩欧美三级三区| 国产精品蜜桃在线观看 | 中国美女看黄片| 中国美白少妇内射xxxbb| 黑人高潮一二区| 晚上一个人看的免费电影| 一进一出抽搐gif免费好疼| 性插视频无遮挡在线免费观看| 亚洲精华国产精华液的使用体验 | 黄色一级大片看看| 亚洲av一区综合| 黄色欧美视频在线观看| 国产精品久久久久久精品电影| 久久精品久久久久久久性| 精品一区二区三区人妻视频| 毛片一级片免费看久久久久| 国产精品乱码一区二三区的特点| 深爱激情五月婷婷| 国产在线男女| 亚洲精品粉嫩美女一区| 老师上课跳d突然被开到最大视频| 91aial.com中文字幕在线观看| 久久久色成人| av.在线天堂| 最近中文字幕高清免费大全6| 欧美最黄视频在线播放免费| 99久久精品热视频| 久久精品久久久久久久性| 日本黄大片高清| 丝袜美腿在线中文| 亚洲最大成人中文| 少妇裸体淫交视频免费看高清| 99热精品在线国产| 久久久久久伊人网av| 免费看av在线观看网站| 亚洲av中文av极速乱| 给我免费播放毛片高清在线观看| h日本视频在线播放| 最近视频中文字幕2019在线8| 国产综合懂色| 日本色播在线视频| 麻豆乱淫一区二区| 91久久精品电影网| 国产片特级美女逼逼视频| av福利片在线观看| 嫩草影院入口| 一个人看的www免费观看视频| 成人漫画全彩无遮挡| 99精品在免费线老司机午夜| 久久中文看片网| 成熟少妇高潮喷水视频| 成人三级黄色视频| 色视频www国产| 一区二区三区高清视频在线| 日本欧美国产在线视频| 精品熟女少妇av免费看| 最后的刺客免费高清国语| 黄色配什么色好看| 岛国在线免费视频观看| 色5月婷婷丁香| 国产精品三级大全| 一级毛片aaaaaa免费看小| 久久中文看片网| 免费黄网站久久成人精品| 国产午夜精品久久久久久一区二区三区| 国产老妇伦熟女老妇高清| 人体艺术视频欧美日本| 一个人观看的视频www高清免费观看| 99久久精品一区二区三区| 久久久久久伊人网av| 欧美另类亚洲清纯唯美| 两个人的视频大全免费| 日日摸夜夜添夜夜爱| 亚洲成人久久爱视频| 日日摸夜夜添夜夜添av毛片| 91久久精品电影网| 男女边吃奶边做爰视频| 人妻系列 视频| www.av在线官网国产| 午夜a级毛片| 亚洲四区av| 国产精品野战在线观看| 人人妻人人澡欧美一区二区| 校园人妻丝袜中文字幕| 女人被狂操c到高潮| 欧美不卡视频在线免费观看| 国产精品免费一区二区三区在线| 亚洲四区av| 亚洲成av人片在线播放无| av免费观看日本| 18禁在线无遮挡免费观看视频| 日本黄色片子视频| 久久鲁丝午夜福利片| 国产成人91sexporn| 99riav亚洲国产免费| 丝袜喷水一区| 国产精品综合久久久久久久免费| 欧美性感艳星| 久久精品国产亚洲av天美| 全区人妻精品视频| 久久婷婷人人爽人人干人人爱| 国产亚洲精品久久久com| 成人永久免费在线观看视频| 18禁黄网站禁片免费观看直播| 91在线精品国自产拍蜜月| 校园春色视频在线观看| 在线观看免费视频日本深夜| 亚洲av.av天堂| 国产片特级美女逼逼视频| 日韩一区二区视频免费看| 边亲边吃奶的免费视频| 乱人视频在线观看| 三级毛片av免费| 成人高潮视频无遮挡免费网站| 男女做爰动态图高潮gif福利片| 在线观看午夜福利视频| 欧美一级a爱片免费观看看| 日韩一本色道免费dvd| 亚洲精品久久国产高清桃花| 性欧美人与动物交配| 亚洲五月天丁香| 三级毛片av免费| av黄色大香蕉| 亚洲四区av| 九色成人免费人妻av| 久久99热这里只有精品18| 国产午夜精品论理片| 一个人观看的视频www高清免费观看| 亚洲成av人片在线播放无| 少妇熟女aⅴ在线视频| 哪个播放器可以免费观看大片| 最近手机中文字幕大全| 老师上课跳d突然被开到最大视频| 网址你懂的国产日韩在线| 五月伊人婷婷丁香| 精品免费久久久久久久清纯| 蜜桃亚洲精品一区二区三区| 欧美日本视频| 只有这里有精品99| 中国国产av一级| 美女 人体艺术 gogo| 看片在线看免费视频| 日本五十路高清| 亚洲欧美日韩高清专用| 国产精品乱码一区二三区的特点| 欧美精品国产亚洲| 久99久视频精品免费| 亚洲在线自拍视频| 12—13女人毛片做爰片一| 国国产精品蜜臀av免费| 亚洲欧美清纯卡通| 中文字幕人妻熟人妻熟丝袜美| 婷婷亚洲欧美| 在线观看一区二区三区| 亚洲精品色激情综合| 国内精品宾馆在线| 国产色爽女视频免费观看| 国产精品一区二区性色av| 亚洲av成人av| 三级毛片av免费| 97超碰精品成人国产| 中文字幕av成人在线电影| 亚洲av免费在线观看| 在线观看午夜福利视频| 99久久人妻综合| 黄色配什么色好看| 国产高清激情床上av| 插阴视频在线观看视频| 男女那种视频在线观看| 1000部很黄的大片| 99在线人妻在线中文字幕| 婷婷色av中文字幕| 在线播放无遮挡| 免费观看a级毛片全部| 18禁裸乳无遮挡免费网站照片| 99热这里只有是精品50| 久久热精品热| 在线观看美女被高潮喷水网站| 人妻久久中文字幕网| av福利片在线观看| 亚洲18禁久久av| 在线免费观看不下载黄p国产| 内射极品少妇av片p| 婷婷色综合大香蕉| 欧美在线一区亚洲| 国产综合懂色| 亚洲精品乱码久久久久久按摩| 欧美性感艳星| 内地一区二区视频在线| 色5月婷婷丁香| 只有这里有精品99| 欧美性猛交╳xxx乱大交人| 男人狂女人下面高潮的视频| 三级经典国产精品| 深夜精品福利| 非洲黑人性xxxx精品又粗又长| 国产伦一二天堂av在线观看| 亚洲欧洲国产日韩| 国产亚洲av嫩草精品影院| 亚洲欧美日韩高清在线视频| 久久久久久九九精品二区国产| 国产精品精品国产色婷婷| 国产美女午夜福利| 亚洲成av人片在线播放无| 久久婷婷人人爽人人干人人爱| 亚洲国产色片| 亚洲精华国产精华液的使用体验 | 国产一区二区亚洲精品在线观看| 非洲黑人性xxxx精品又粗又长| 国产午夜精品论理片| 国产蜜桃级精品一区二区三区| 亚洲欧美中文字幕日韩二区| 麻豆一二三区av精品| 99久久精品热视频| 国产精品久久久久久精品电影小说 | 看免费成人av毛片| 国产精品国产三级国产av玫瑰| 国产美女午夜福利| 两性午夜刺激爽爽歪歪视频在线观看| 九九久久精品国产亚洲av麻豆| 亚洲欧美清纯卡通| 国产成人午夜福利电影在线观看| 亚洲成人中文字幕在线播放| 午夜老司机福利剧场| 免费av观看视频| 26uuu在线亚洲综合色| 亚洲av免费高清在线观看| av在线播放精品| 亚洲天堂国产精品一区在线| 国产高潮美女av| 国产黄片视频在线免费观看| 久久久久久久亚洲中文字幕| 久久久午夜欧美精品| 亚洲精品影视一区二区三区av| 大又大粗又爽又黄少妇毛片口| 久久精品影院6| 观看美女的网站| 国产一级毛片在线| 两性午夜刺激爽爽歪歪视频在线观看| 久久中文看片网| 国产亚洲欧美98| 久久精品91蜜桃| 亚洲av第一区精品v没综合| 国产成人精品一,二区 | 非洲黑人性xxxx精品又粗又长| 夫妻性生交免费视频一级片| 国产精品日韩av在线免费观看| 精品人妻一区二区三区麻豆| 久久99精品国语久久久| 给我免费播放毛片高清在线观看| 亚洲最大成人中文| 欧美性猛交黑人性爽| 老女人水多毛片| 国产不卡一卡二| 国产日本99.免费观看| 亚洲一区二区三区色噜噜| 一个人观看的视频www高清免费观看| 亚洲成av人片在线播放无| 五月玫瑰六月丁香| 色尼玛亚洲综合影院| 久久99蜜桃精品久久| 麻豆国产av国片精品| 成人亚洲欧美一区二区av| 午夜亚洲福利在线播放| kizo精华| 一级二级三级毛片免费看| 精品人妻一区二区三区麻豆| 国产真实伦视频高清在线观看| 亚洲欧美日韩高清专用| 国产黄片美女视频| 天美传媒精品一区二区| 少妇的逼水好多| 久久人人爽人人爽人人片va| 中文字幕免费在线视频6| 国产av在哪里看| 在线免费十八禁| 插逼视频在线观看| 国产高清有码在线观看视频| 亚洲人成网站在线观看播放| 精品午夜福利在线看| 高清在线视频一区二区三区 | 久久久久久久久久成人| 一级毛片电影观看 | 插阴视频在线观看视频| 老熟妇乱子伦视频在线观看| 欧美性猛交╳xxx乱大交人| 亚洲三级黄色毛片| 丰满的人妻完整版| 国产精品蜜桃在线观看 | 亚洲自偷自拍三级| 国产av一区在线观看免费| 亚洲成人中文字幕在线播放| 午夜激情福利司机影院| 国产极品精品免费视频能看的| 日韩高清综合在线| 99久久九九国产精品国产免费| 久久综合国产亚洲精品| 亚洲成av人片在线播放无| 亚洲内射少妇av| 亚洲在线观看片| 免费电影在线观看免费观看| 亚洲国产精品国产精品| 日韩精品有码人妻一区| 久久草成人影院| 久久久午夜欧美精品| 国产美女午夜福利| 国产女主播在线喷水免费视频网站 | 欧美3d第一页| 2022亚洲国产成人精品| 51国产日韩欧美| 国产爱豆传媒在线观看| 精品欧美国产一区二区三| 在线天堂最新版资源| 欧美xxxx黑人xx丫x性爽| 三级经典国产精品| 嫩草影院新地址| 亚洲精华国产精华液的使用体验 | 国产又黄又爽又无遮挡在线| 99久久精品热视频| 欧美另类亚洲清纯唯美| 欧美日韩精品成人综合77777| 青春草国产在线视频 | 亚洲精品色激情综合| 老女人水多毛片| 亚洲av男天堂| 精品人妻一区二区三区麻豆| 久久中文看片网| 国内久久婷婷六月综合欲色啪| 国内揄拍国产精品人妻在线| 亚洲av电影不卡..在线观看| av免费观看日本| 免费看日本二区| 2021天堂中文幕一二区在线观| 少妇人妻精品综合一区二区 | 久久亚洲精品不卡| 亚洲久久久久久中文字幕| 18禁黄网站禁片免费观看直播| 如何舔出高潮| av国产免费在线观看| 美女高潮的动态| 2021天堂中文幕一二区在线观| 99热全是精品| 久久久午夜欧美精品| 少妇人妻精品综合一区二区 | 亚洲自偷自拍三级| 亚洲人成网站高清观看| 97超视频在线观看视频| 久久久久久久亚洲中文字幕| 日韩三级伦理在线观看| 国产成人精品一,二区 | 精品久久久久久久久久久久久| 菩萨蛮人人尽说江南好唐韦庄 | 一夜夜www| 日韩,欧美,国产一区二区三区 | 久久久精品94久久精品| 久久这里有精品视频免费| 久久久久久国产a免费观看| 麻豆一二三区av精品| 久久婷婷人人爽人人干人人爱| 蜜桃久久精品国产亚洲av| 久久精品国产亚洲av香蕉五月| 久久精品国产亚洲av天美| 99久久精品热视频| 联通29元200g的流量卡| 国产成人精品婷婷| 午夜视频国产福利| 偷拍熟女少妇极品色| 国产成人freesex在线| 观看美女的网站| av福利片在线观看| 中国美女看黄片| 国产精品久久久久久久电影| 国产精华一区二区三区| 亚洲第一区二区三区不卡| 久久久国产成人免费| 在线观看av片永久免费下载| 一级黄片播放器| 日韩欧美精品v在线| 欧美另类亚洲清纯唯美| 久久久久久久久大av|