999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

含天然氣水合物土微觀力學特性研究進展1)

2021-11-10 03:43:22趙亞鵬劉樂樂劉昌嶺吳能友
力學學報 2021年8期
關鍵詞:實驗模型研究

趙亞鵬 劉樂樂 孔 亮?,3) 劉昌嶺 吳能友

* (青島理工大學土木工程學院,山東青島 266033)

? (中國地質調查局青島海洋地質研究所,自然資源部天然氣水合物重點實驗室,山東青島 266071)

** (青島海洋科學與技術國家實驗室海洋礦產資源評價與探測技術功能實驗室,山東青島 266071)

?? (青島理工大學理學院,山東青島 266033)

引 言

天然氣水合物(簡稱水合物)是由天然氣分子與水分子在高壓、低溫環境下形成的非化學計量籠型結晶物質[1],廣泛分布于全球近海大陸架及陸地永久凍土地帶[2],資源儲量巨大,被認為是21 世紀最具潛力的替代能源之一[3].自20 世紀開始,世界各國陸續成立了專門的水合物研究機構并開展了相關研究[4],在此基礎上,前蘇聯、加拿大、美國及日本等國家[5-8]相繼進行了水合物試采,證明了水合物開采的技術可行性.與此同時,經過20 余年的不懈努力,我國在南海海域發現了儲量高達800 億噸油當量的水合物資源[9],并于2017 和2020年在南海北部海域兩次成功試采水合物[10-11],創造了多項世界記錄,率先實現了泥質粉砂型水合物資源的安全可控開采.

水合物資源的開發前景主要受開采技術條件及開發成本的控制,應以實際工程地質條件為基礎,結合儲層特性選取合理的開采方式[12].當前水合物開采方法主要有熱激、降壓、注入化學試劑、二氧化碳置換以及多方法聯合等[13-17],雖然開采方法各異,但是均需原位分解水合物[18-19].水合物分解導致沉積物顆粒膠結程度降低,并可能產生超靜孔壓,儲層力學強度隨之弱化,進而引發各種工程地質災害,如海底滑坡、海底濁流、井壁失穩及平臺傾覆等[20-23].因此,厘清開采過程中含水合物土的力學特性演化規律是實現工程地質風險可控的重要理論基礎,是實現水合物資源安全、高效、持續開采的前提條件之一[24-25].

國內外學者圍繞含水合物土的力學特性開展了大量的研究工作,研究態勢呈現出多尺度融合、多方法聯用、多學科交叉等特征.目前實驗研究主要以高壓低溫水合物三軸剪切實驗為主[26],輔以電學、光學、聲學等測量手段[27-29],再結合理論分析、數值模擬[30-31]等研究方法,系統研究含水合物土的力學特性演化機制.研究表明,含水合物土的力學特性由其內部結構控制,特別是水合物的微觀行為扮演了極其重要的角色,越來越多的學者著眼于水合物微觀力學領域[32],試圖從微觀結構演化解釋宏觀力學現象,“宏微觀相結合”的科學范式成為當前乃至未來含水合物土力學特性研究的主流[33].因此,本文從水合物晶體(分子尺度)、水合物與土顆粒界面(顆粒尺度)、含水合物土(巖心尺度)3 個尺度對含水合物土微觀力學特性研究現狀進行綜述,歸納微觀力學特性研究的最新進展,闡述主要的微觀測試技術,梳理出當前含水合物土微觀力學研究的不足并給出針對性建議,以期為含水合物土力學特性研究提供新的思路與方向.

1 水合物晶體力學性質研究

1.1 晶體結構類型

水合物晶體的“主體”水分子之間由氫鍵連接并形成“籠型”結構,而“客體”天然氣分子則通過范德華力填充于“籠”中[34].學者們現已發現多達100 余種水合物,按結構將其分為3 類,即sI,sII 及sH 型[35].不同類型水合物具有不同的晶體結構,在抗壓強度和體積模量等力學性質上也表出明顯的差異[36-37].sI,sII 和sH 型水合物所能容納的烴類氣體分子依次增大,而最終形成何種類型水合物則與氣體組分及其含量有關.不同結構類型水合物對比可見表1.

表1 不同類型天然氣水合物晶體結構對比[36,38]Table 1 Comparison of different structures in natural gas hydrate crystals[36,38]

1.2 研究方法

1.2.1 實驗研究

考慮到晶粒尺寸對多晶固體力學性能的顯著影響[39],學者們預猜想多晶水合物在力學特性上也將受到晶粒大小的影響[40].針對該問題,早期研究對比分析了海洋水合物與合成水合物的晶粒尺寸分布[41],并進行了毫米尺度下二氧化碳(CO2)水合物、甲烷(CH4)水合物以及四氫呋喃(C4H8O)水合物的抗拉強度研究[42],印證了晶粒尺寸確實會對水合物晶體力學特性尤其是塑性應變產生重要的影響.Atig 等[43]介紹了一種非常規的水合物微觀力學特性研究方法,探討了多晶甲烷水合物在微米尺度下的拉伸特性.該研究首先在薄玻璃毛細管中生成水合物,然后將非接觸的熱應力施加于水合物外殼(圖1),進而研究了微米級別的多晶水合物抗拉特性.研究結果表明,溫度和退火時間可控制水合物晶粒大小,水合物脆性隨退火時間的延長而增強,而極限強度及應變則顯著降低;顆粒結構的冰晶體并不是多晶甲烷水合物的理想相似物,海域甲烷水合物的抗拉強度可能遠低于目前業界的認識.

圖1 非接觸應力施加步驟[43]Fig.1 Principal application steps of contactless stress[43]

水合物顆粒間的黏聚力與周圍水的接觸時間、接觸力及界面能成正比,并隨溫度的升高而增大[44],水合物顆粒之間的液橋力是引起顆粒黏附力的主要原因[45].當含水合物土受到外力作用而發生顆粒旋轉、滾動時,首先需要克服顆粒液橋力,這一點在水合物力學特性的原子力顯微鏡(atomic force microscopy,AFM)研究中已得到證明.液橋力示意圖及AFM 實驗液橋模型如圖2 所示.

圖2 液橋力示意圖及AFM 實驗液橋模型[46-47]Fig.2 Schematic of liquid bridge force and the model of liquid bridge for AFM experiment[46-47]

根據相關文獻,即使在低于相平衡的溫度下,水合物表面仍然存在納米級別的似液層[48]以及壓力作用下的分解液[49],氣體水合物與砂顆粒之間亦會存在[50].顆粒間不同作用力(范德華力、液橋力、流體作用力等)的理論計算[51],則從定量的角度上進一步證明了液橋力是導致水合物顆粒黏附力的主要原因.而顆粒間液橋力受多因素的影響,以溫度及液橋體積影響為切入點,對環戊烷水合物顆粒間的相互作用研究表明[52],未轉化的水滴對顆粒團簇具有重要影響,液橋的形成能顯著增強顆粒間黏附力,而黏附力、接觸面積則與液滴體積沒有直接的關系;同時,水合物顆粒呈現出與冰顆粒相似的力學行為[53].借助于水合物AFM,學者們用液橋力的公式進行了計算[46],發現探針與水合物間黏附力隨溫度的降低而減小.此外,受環境因素的影響,水合物顆粒與液體水滴呈動態作用關系,兩者不斷轉化的同時伴隨體積的持續變化,這將對顆粒間的液橋力產生顯著影響[54].Aman 等[45,47,55]開展了較為系統的水合物液橋力研究,通過改進的微機械力學裝置提出了一種直接測量顆粒間黏附力的新方法,對氣?水兩相黏附力進行對比分析,進一步探討了不同添加劑對水合物顆粒之間液橋力的影響.

現階段顆粒液橋研究主要集中于原子力顯微鏡等微觀裝置的實驗分析,在進一步擴大研究尺度以及宏觀聯系方面范圍有限,應逐步開展顆粒液橋的宏觀力學影響研究.另一方面,雖然目前已經建立起了一定數量的液橋力微觀理論模型,但是考慮多因素影響的理論模型仍相對匱乏,也有待進一步地深入研究.

1.2.2 分子動力學模擬

在晶體力學性質研究方面,除實驗研究外,分子動力學(molecular dynamics,MD)模擬也是較為常用的研究方法.MD 以牛頓力學為基礎,結合疊加原理,能夠對分子尺度的微觀力學進行模擬研究.

MD 研究表明水合物晶體強度、彈性模量及變形等與客體分子的大小、形狀、極性等密切相關[56],同時還受到應變速率、溫度、客體分子占有率的影響[57].水合物晶體與冰晶體的力學特性相近,一般表現為脆性破壞,而水合物斷裂取決于客體分子類型,變形則受氫鍵的角位移控制.主客體分子關系的MD 研究發現[58],甲烷分子在“大籠”及“小籠”中均呈現各向同性,二氧化碳分子則在“大籠”中表現出明顯的各向異性;而在混合氣體中,無論是客體分子之間,還是主客體分子之間,甲烷的存在都會導致二氧化碳分子分布傾向性的改變.與純水合物相比,冰含量對于含冰水合物強度影響較大,水合物?水合物雙晶結構抗剪強度遠大于水合物?冰雙晶抗剪強度[59].在單晶及多晶水合物力學特性方面,單晶水合物呈現脆性破壞,多晶則為延性破壞,多晶甲烷水合物的力學穩定性與顆粒大小及形態密切相關,這種依賴關系歸因于晶粒尺寸對晶界變形的影響[60].此外,學者們圍繞天然氣水合物晶體,通過MD 方法開展了包括不同氣體分子置換過程[61]、水合物生長成核[62-63]、分解動力學[64]及分子作用模型[65]等方面的研究,進一步補充和豐富了水合物晶體分子尺度研究.水合物的典型MD 模擬如圖3 所示.

圖3 甲烷水合物分子結構[59]Fig.3 Molecular structure of methane hydrate[59]

相對而言,MD 模擬是一種非常規的研究手段,可以從分子及原子尺度進行微觀、介觀到宏觀的物理力學性質解釋.在水合物微觀力學特性研究尤其是在水合物成核、分解及熱力學方面,起到了極大地推動作用.而當前水合物MD 模擬研究主要基于均勻成核,而真實自然環境中存在非均勻成核情況,如何合理模擬非均勻成核帶來的影響將是MD 研究的一大挑戰.近期水合物晶體MD 模擬研究匯總可見表2.

表2 近期水合物晶體MD 模擬研究匯總Table 2 Summary of recent MD research on hydrate crystals

2 水合物與土顆粒界面力學性質研究

將水合物晶體結構尺度進一步擴大,晶體團簇形成顆粒,水合物顆粒與沉積物顆粒之間按照接觸形態的不同可分為4 種主要賦存模式[67](如圖4 所示),圖4 中GH 為水合物,SP 為土顆粒.圖4(a)為孔隙填充模式.水合物在沉積物顆粒孔隙之間生成,與沉積物顆粒無實質接觸[68],此時水合物飽和度較低;圖4(b)為顆粒支撐模式.在顆粒填充模式基礎上,水合物顆粒進一步擴大并與沉積物顆粒接觸,共同承擔應力傳遞,此時水合物飽和度達到25%~40%[69];圖4(c)為膠結模式.與顆粒填充模式不同,膠結模式水合物首先在沉積物顆粒接觸處生成,并將周圍沉積物顆粒“黏結”在一起,形成整體骨架[70];圖4(d)為包裹模式.該模式下水合物將沉積物顆粒完全包裹.

圖4 粗粒土中水合物賦存模式[67]Fig.4 Idealized gas hydrate morphologies in coarse-grained soils[67]

以上4 種模式主要存在于粗粒沉積物中,除孔隙填充模式外,其余3 種模式均會對力學特性產生顯著影響[71].在“限氣”或“富水”情況下易形成填充模式,在“限水”或“富氣”情況下易形成包裹及膠結模式[72-73].細粒沉積物中則以結核狀、脈狀、塊狀、透鏡狀等模式存在[74-75](見圖5),生成過程傾向于“擠開”沉積物發生“顆粒替代”,這與粗粒沉積物區別明顯.上述賦存模式可進一步分為粗粒沉積物中均勻分布的孔隙填充型水合物以及細粒沉積物中非均勻分布的裂隙填充型水合物[76-77].對于含不同賦存類型水合物土的力學性質,很大程度上受到水合物與土顆粒界面力學行為的影響與控制.

圖5 含水合物細粒土實物巖心[77]Fig.5 Hydrate-bearing fine-grained soils in nature[77]

2.1 研究方法

2.1.1 模擬實驗研究

水合物與土顆粒界面力學性質的實驗研究以計算機斷層掃描(computed tomography,CT)聯用的三軸剪切實驗為代表,還包括部分基于微觀實驗裝置、微觀測試技術的顆粒界面關系研究.Lei 等[78]的CT 聯用三軸實驗表明,水合物主要通過限制顆粒旋轉強化顆粒骨架,進而實現沉積物強度及剛度的增大,而應變軟化的主要原因在于水合物與砂顆粒之間裂紋的產生;剪切應力會引起顆粒破碎并產生破壞裂縫,水合物則“流入”剪切裂縫或被尖銳的破碎顆粒所刺破.針對剪切過程的局部變形問題[79-80],借助于CT 圖像,學者們確定了試樣離散位置的局部剛度,并發現局部剛度與水合物飽和度成正相關關系,其中不含水合物的沉積物樣品表現出延性破壞,含水合物樣品則呈現出脆性破壞.基于微米尺度圖像的局部變形(剪切面、顆粒變形及旋轉等)量化分析[81],發現剪切帶內沉積物的結構變化是基質顆粒和水合物顆粒旋轉運動的結果,剪切帶厚度隨水合物飽和度的增大而減小.此外,水合物首先在顆粒表面發生分解,隨后向顆粒孔隙發生轉移;而水合物分解會引起沉積物組構的變化,進一步增強砂土顆粒分布的各向異性,并造成顆粒支撐及膠結效應的顯著降低,顆粒間的均勻分布則可顯著增強砂粒體系的力鏈穩定性[82].表3 給出了近期水合物CT 三軸實驗匯總.

表3 近期含水合物土CT 聯用三軸剪切實驗匯總Table 3 Summary of current triaxial shearing tests combined with CT conducted on hydrate-bearing soils

除CT 三軸實驗外,基于微觀力學測量裝置的直接測量結果表明[44,86],界面黏附力與顆粒接觸時間及載荷大小有關,在接觸時間一定的條件下,黏附力隨載荷的增大而增大;顆粒摩擦特性方面,水合物顆粒與沉積物顆粒之間摩擦系數在不同載荷下呈現出不同的變化規律,當載荷較小時,隨載荷增大,水合物顆粒與沉積物顆粒間摩擦系數亦不斷增大.根據水合物的AFM 研究[87],顆粒間摩擦特性的差異在于顆粒表面性質或顆粒接觸關系的不同,如顆粒表面粗糙度、顆粒尺寸、顆粒間似液層厚度等;水合物的不定形特征主要受低溫及小尺寸晶粒的誘導,水合物表面形貌及粗糙度受環境溫度、接觸介質屬性及晶粒尺寸所控制.不同個體分子的水合物抗拉實驗則進一步表明[88],水合物在基質上成核以及水合物與基質整體接觸時,表面粗糙程度對宏觀力學特性的影響是不同的.此外,與水合物晶體顆粒類似,水合物與沉積物顆粒間黏聚力也會受到液橋力的影響,且沉積物基質表面對液橋力具有積極影響[45].

2.1.2 數值模擬研究

與有限元解決連續問題不同,離散元法(discrete element method,DEM)的顆粒流程序(particle flow code,PFC)不受變形量限制,能夠對剪切過程中的開裂、分離等現象進行模擬,特別適用于解決非連續介質力學問題[89],被廣泛應用于微觀力學分析,并被認為是聯系宏?微觀力學的有效橋梁[32].

根據PFC 研究,顆粒界面間的黏結作用是含水合物土宏觀力學響應的微觀本質,這種黏結作用會對受力過程中的顆粒滾動產生顯著影響,限制顆粒翻滾,進而影響宏觀性質.CT 與掃描電子顯微鏡(scanning electron microscope,SEM)研究[90-91]表明,沉積物顆粒之間具有不同的顆粒形狀,呈現不規則特性,顆粒形狀對土體的影響則主要體現在法向接觸力偏離形心引起的滾動力矩作用[92].同時不同形狀顆粒具有不同的抗滾動系數,基于顆粒抗滾動作用的微觀數值模擬表明[93],抗滾動系數的增大雖然會導致配位數的降低、孔隙率的增大,但是由于平均接觸力的增大,顆粒之間將會形成更強的力鏈,進而引起宏觀力學特性的增強.目前室內實驗大多采用氣飽和法合成水合物,生成的水合物以顆粒包裹型為主,相關的PFC 數值模擬研究表明,從顆粒界面關系的角度出發,包裹型水合物存在3 種微觀作用機制[94],即顆粒黏結、顆粒擴大以及顆粒成角(黏結斷裂時,存在粗糙斷裂面,增大顆粒粗糙度).基于3 種微觀機制的真三軸DEM 模擬[95]表明,包裹型水合物顆粒間的黏結作用相對較弱,真三軸剪切過程的力學特性接近于無黏結顆粒松散材料.此外,水合物剪切過程中在剪切帶內外存在明顯的顆粒演化差異[96],剪切帶內部黏結大量破壞,顆粒發生轉動,局部孔隙率不斷增大;剪切帶外部黏結破壞則相對較少,顆粒幾乎不發生轉動.這也進一步表明水合物與沉積物顆粒界面間的黏結關系是水合物存在而導致的微觀效應,其會造成沉積物剪切破斷過程中的微觀顆粒旋轉抑制,并呈現區域差異性,最終導致沉積物宏觀力學特性的演變.

PFC 對水合物力學特性的有效模擬依賴于其內部的計算模型,而水合物顆粒與基質顆粒之間的界面作用關系,尤其是黏結效應則是PFC 計算模型建立的基礎.不同類型水合物與沉積物顆粒之間具有不同的黏結厚度,這種界面關系將會呈現出宏觀上的差異性特征[97].針對該問題,借助于SEM 圖片,可獲得水合物飽和度與黏結厚度之間的函數關系,進而建立考慮水合物黏結厚度的微觀力學膠結模型[98].在多場耦合研究方面,基于實驗結果,還可建立考慮應力、溫度影響的黏結接觸離散元模型[99],能夠有效反映應變局部化特征;以及考慮水合物黏結厚度的熱與流體力學耦合模型[100],對該模型的雙軸壓縮實驗表明,隨著顆粒間黏結厚度的增大,水合物峰值強度及剪脹率均顯著增大.將計算流體動力學與DEM 相結合,學者們提出了一種考慮水合物速率依賴性的黏結接觸模型[101](如圖6 所示),基于該模型的不排水循環剪切實驗表明,力鏈、接觸組構以及平均旋轉速率等微觀變量與宏觀力學特性之間具有緊密的聯系.此外,部分學者認為水合物的存在會導致摩擦性的增強,而不一定與黏結有關,據此提出了含水合物土的非黏結微觀模型[102].模型中將水合物表示為固體顆粒,即使在小應變下也會產生骨架效應.顆粒粒徑與顆粒摩擦力之間存在著平衡關系,利用與孔徑分布相關的附加因子可以實現兩者的統一.

2.2 研究結果

基于實驗結果,學者們試圖從水合物與沉積物顆粒間界面關系,尤其是其在剪切過程中的動態演化規律對宏觀力學特性進行解釋.現有實驗數據表明,水合物的存在會對沉積物硬化、軟化、剪縮、剪脹等特性產生顯著影響,顆粒尺度分析認為這是水合物顆粒、沉積物顆粒及孔隙空間相對關系的宏觀體現[103-104].

水合物的生成會占據孔隙空間,這種“填充”作用導致沉積物整體骨架致密性的提高,同時水合物對砂顆粒具有顯著的“黏結”作用,兩者共同對力學特性產生影響.而不同賦存模式下含水合物土的“填充”與“黏結”作用也會存在明顯差異,其受水合物飽和度的控制,內在本質則是微觀組構關系.針對剪切過程中的含水合物土,隨著剪切應力的升高,顆粒間接觸增多,并不斷摩擦、擠壓及旋轉,最終發生顆粒破碎以及“黏結”破壞.即剪切過程的實質是“填充”與“黏結”作用控制下的顆粒關系重新排列.對于密實沉積物而言,往往發生剪脹,松散沉積物則對應剪縮[105].

剪切過程的顆粒尺度作用機理可見圖7.當水合物飽和度較低時,水合物對沉積物的強度影響較小,此時水合物與砂顆粒之間接觸較少,“黏結”作用不明顯[106-107],而水合物自身的“膠結”強度很容易被剪切應力所突破[108].隨著水合物飽和度的增大,水合物與砂顆粒發生接觸,“黏結”作用增強,水合物成為顆粒骨架并與砂顆粒共同承擔應力傳遞.此時水合物的存在會顯著限制基質顆粒在剪切應力下的滾動及重新排列.在此基礎上,水合物強度演化受水合物“膠結”強度、砂顆粒與水合物之間“黏結”強度相對大小的控制.當“黏結”強度大于“膠結”強度時,首先在水合物中產生剪切面,沉積物強度演化由水合物力學特性控制,這種情況很可能發生在表面粗糙或具有較低比表面積沉積物中;當“膠結”強度大于“黏結”強度時,水合物與砂顆粒界面處將產生剪切面,并在剪切面形成大量水合物顆粒,此種情況很可能發生在表面光滑或具有較高比表面積沉積物中[69].

圖7 顆粒尺度下含水合物土剪切機理[105]Fig.7 Particle-scale shearing mechanism of hydrate-bearing soils[105]

針對水合物顆粒的剝離與破碎問題,一些學者認為含水合物土在拉伸及剪切過程中分別具有不同的破壞模式[109].一種假想的微觀機制如圖8 所示,當受拉應力作用時,顆粒剝離或拉伸破壞主要取決于界面性質,并呈現出兩種破壞模式,分別為沉積物顆粒與水合物顆粒之間的剝離以及水合物自身內部的斷裂.失效的水合物與固體顆粒之間沒有進一步的相互作用,兩種模式如圖8(a)所示.而在剪切破壞中,破碎顆粒將繼續與其他顆粒產生相互作用,宏觀上表現為峰后強度[69,110-111].與拉伸破壞相比,剪切破壞存在混合機制,即顆粒剝離過程中伴隨水合物的斷裂,相比純顆粒剝離而言其體變更小.

圖8 含水合物土拉伸及剪切破壞機制(改自文獻[109])Fig.8 Failure modes of tension and shearing in hydrate-bearing soils (modified from Ref.[109])

近期,Wu 等[85]通過CT 三軸實驗對水合物的拉伸及剪切過程進行了更深一步的探討分析.CT 圖像(圖9)表明,無論是哪種類型的水合物,在剪切過程中,水合物顆粒均呈現整體移動,顆粒間相對位置關系幾乎不變,而不含水合物的砂顆粒即使在較小軸向應變下也呈現無序移動.這無疑進一步印證了水合物對沉積物整體強度的顯著增強作用.并且水合物對顆粒重新排列的抑制作用在應力應變的線性區域(εa≤2%)效果明顯;隨著軸向應變增大,應力應變曲線進入塑性及屈服階段,水合物脫落、顆粒破碎,并出現了明顯的剪切帶,最終在剪切作用下發生了顆粒的重新排列.

圖9 含水合物土剪切過程中截面CT 圖像[85]Fig.9 CT images of hydrate-bearing soils during shearing[85]

3 含水合物土力學性質研究

3.1 離散元數值模擬

PFC 對水合物的有效模擬以制樣方法為基礎,重點在于突出水合物及沉積物顆粒之間的相互作用關系,進而獲取這種微觀關系所導致的宏觀力學特性.就黏結效應而言,現有研究發現微觀黏結參數對宏觀力學特性的影響存在差異[112],隨黏結強度的增大,峰值強度及彈性模量均增大;隨黏結半徑的增大,彈性模量增大而峰值強度基本不變.在抗滾動方面,由于黏結存在而導致的抗滾動作用會顯著提高含水合物土的抗剪強度以及剪脹性,但是其存在一個閾值,當超過閾值后,這種增強作用將逐漸減小[93];并且不同的抗滾動作用會形成不同的微觀力鏈分布(如圖10 所示),這也將對沉積物宏觀力學特性產生顯著影響.含水合物土的PFC2D 數值模擬表明[113],剪切帶的形成與黏結的大量破壞有關,并與應變局部化、微觀變量局部化密切相關.以黏土質沉積物為研究對象,并考慮顆粒間毛細水效應的模擬研究發現[114],含水合物土具有剪縮及應變硬化的力學特性;隨水合物飽和度的增大,顆粒間毛細水合物將產生一種鍵合效應,其隨顆粒間隙的增大而線性減小.根據多類型實驗(三軸壓縮、各向同性壓縮、恒應力比壓縮及真三軸壓縮)的水合物宏微觀力學特性研究[115],顆粒包裹型水合物力學特性表現為弱黏結砂,顆粒膠結在加載初期就破裂,而在加載后基質砂顆粒的摩擦特性隨顆粒形貌的變化而變化.對循環載荷下的含水合物土DEM 模擬表明[116],沉積物的破壞模式取決于循環載荷幅值以及水合物飽和度,組構的各向異性存在閾值,超過閾值沉積物結構則會趨于破壞.此外,針對含水合物土的黏聚力及內摩擦角力學參數,PFC 模擬研究發現[117],黏聚力和內摩擦角隨剪脹量的增大而減小,應變局部化與黏結破壞、接觸力鏈、顆粒速度等微觀局部變量密切相關.當水合物飽和度增大時,水合物抗剪強度、黏聚力及體縮均明顯增大,而摩擦角基本不變[96,100,118],這與宏觀三軸實驗結果一致[119].

圖10 不同抗滾動系數下含水合物土內力鏈分布[93]Fig.10 Force chain distribution in hydrate-bearing soils under different rolling resistance coefficients[93]

微觀參數對PFC 模擬結果至關重要,但是由于離散元顆粒與真實顆粒之間的差異,微觀參數的確定較為困難,且不能直接獲得.目前主要通過“試算”的方式確定微觀參數[89,112],首先以宏觀水合物三軸實驗為依據,在一定范圍內賦予接觸模型微觀力學參數,根據模擬結果不斷調試,最終獲得有效反應宏觀力學特性的微觀力學參數[120],據此開展一系列的微觀數值模擬研究.同時基于前述的PFC 內部計算模型分析,可知,構建有效反應水合物力學特性的內部計算模型是水合物PFC 研究的另一難點與基礎.針對該問題,蔣明鏡等[121-122]開展了早期的探索研究,并進行了有關模型建立與參數反演的深入分析,研究成果對水合物PFC 研究起到了很好的推動作用.

DEM 以PFC 為代表在巖石力學及巖土工程領域得到了廣泛應用,當前水合物力學研究所涉及的諸多力學問題都可劃歸為巖土力學問題[33].同時PFC 可實現水合物從微米顆粒尺度[123]到海底滑坡千米尺度[124]的有效模擬,因此可以預見,未來在水合物宏?微觀相結合研究尤其是微觀力學特性方面,PFC 還將繼續發揮出更大的作用.

3.2 本構模型研究

本構模型研究主要以力學實驗結果為依據,結合彈塑性本構建模理論及相關假設進行模型建立.在水合物方面,目前較多的是基于水合物存在對沉積物宏觀力學特性的影響而開展相關研究,例如通過含水合物土的三軸剪切實驗,對不同水合物飽和度的應力應變關系進行分析,進而建立起適用于中細粒含水合物土的本構模型[125].微觀力學本構則主要圍繞水合物與沉積物之間的微觀作用關系、力學作用機理等進行模型構建,如考慮水合物賦存模式、填充及黏結效應、基質顆粒作用關系、顆粒破碎及損傷等.

水合物的微觀賦存模式對含水合物土的宏觀力學特性具有顯著影響,這已經成為普遍共識.基于這一考慮角度,通過引入水合物“有效飽和度”的概念,并重新定義有效應力,可以建立起考慮水合物微觀賦存模式的彈塑性本構模型[126],能夠有效反映水合物微觀賦存模式及水合物飽和度對含水合物土宏觀力學特性的影響,對剛度、強度、應變軟化等力學特性的演變均能較好地描述[127].針對水合物的微觀黏結機理,在經典的臨界狀態模型框架內,通過引入非線性黏結與線性脫結規律,并將臨界狀態線表示為水合物飽和度的函數,建立起考慮水合物微觀黏結機理以及黏結強度的本構模型[128],能夠充分反映水合物飽和度、圍壓以及密度對含水合物土力學特性的影響;此外,引入次加載概念,用增強因子表征水合物的微觀膠結效應,同時結合顆粒損傷,則可建立起能夠有效預測不同孔隙習性(如孔隙填充、孔隙膠結等)水合物力學行為的新型本構模型[129].水合物的存在會造成含水合物土微觀孔隙結構的改變,部分學者將這種孔隙改變歸結為填充及膠結效應,通過引入壓硬性參量,建立起考慮水合物填充及膠結效應的彈塑性本構模型[130].沉積物基質與水合物顆粒之間存在明顯的界面關系,這將造成宏觀受力過程中,微觀力學響應的差異,即兩者具有不同的力學貢獻,據此,將宏觀受力進行微觀劃分,同時引入“分區應力”的概念,則可建立起能夠有效描述水合物加載及分解過程中的本構模型[131].Fang 等[132]根據塑性滑移理論,將宏觀受力分解為體積響應以及空間分布的一系列微剪切響應(如圖11 所示),進而得到了本構方程,該模型可以綜合描述含水合物土的固結、硬化、軟化、剪脹及非共軸等特性.

圖11 含水合物土的多剪切模型概念框架[132]Fig.11 Illustration of conceptual framework for multishear model of hydrate-bearing soils[132]

針對水合物開采出砂問題,通過不同方法對微觀顆粒運移的本質進行量化研究,發現水合物的存在減緩了應力松弛幅度,據此提出了考慮應力松弛特性的本構模型[133],該模型可實現水合物儲層產氣行為的有效模擬.部分學者[134]認為含水合物土缺乏真正的黏聚力,其力學行為受孔隙尺度的顆粒運動所控制.據此建立起了一種全新的水合物彈塑性本構模型,該模型并不依賴于水合物晶體與沉積物顆粒間的物理黏結,而是依賴于水合物侵入孔隙并對沉積物產生的致密化效應.根據該模型,水合物形成過程中沉積物微觀孔隙體積的減小使其結構變硬,并隨著沉積物密度的增大而產生近似的力學效應.此外一些學者們從損傷力學的角度出發,認為應力?應變關系的塑性階段是內部顆粒發生損傷并逐漸累積的過程[135],損傷過程則可分為微觀結構的快速損傷與完全損傷階段[136],同時結合微元強度、細觀力學機制等微觀假設,建立起了一系列的基于損傷理論的含水合物土本構模型[137-139].其建模思路及建模過程對水合物微觀力學本構模型的建立起到了很好的借鑒作用.

水合物本構模型及微觀力學本構模擬方面已經積累了大量的研究成果,而其主要問題在于模型適用性及參數確定方面,常用的模型參數確定方法包括實驗獲取、擬合分析以及合理假設等[127,130,140-142].絕大多數模型基于特定實驗條件或特定假設而建立,適用范圍有限.在模型參數方面,部分模型的參數眾多,逐一確定參數取值十分困難,很多模型參數并不具有實際物理意義,這無疑進一步增大了確定參數的難度.從建模角度出發,上述本構模型可進一步劃分為3 類,分別為彈塑性模型、統計損傷模型以及其他模型.而宏觀水合物本構模型還包括非線性彈性模型,然而由于該類模型采用彈性框架,在水合物微觀力學特性建模方面受到了較大局限,報道較少.含水合物土微觀本構模型分類及主要特點如表4 所示.

表4 含水合物土微觀本構模型分類及特點Table 4 Classification and characteristics of microscopic constitutive models of hydrate-bearing soils

4 微觀力學實驗測試技術

當前主要通過三軸剪切及直剪實驗對含水合物土力學特性進行研究,尤其在三軸剪切實驗方面成果頗豐[125,143-145].微觀力學實驗方面則可根據其研究方式的不同分為兩類,一類是將宏觀力學實驗與微觀測試技術相結合,通過水合物微觀結構的演變對宏觀力學特性進行分析[106],屬于間接方式;另一類是基于三軸剪切等力學實驗原理,研發專用于微觀力學特性研究的實驗裝置,直接測量水合物及其界面力學參數的測試技術,如原子力顯微鏡、微三軸儀等[46,81].

4.1 間接測試技術

間接測試技術一般只測量研究對象的結構,力學性質需要結合三軸實驗、直剪實驗及聲學實驗等實驗手段.隨著現代測試及分析技術的發展,眾多具有交叉學科屬性的測試技術被廣泛應用于水合物研究,而在微觀力學領域,則以光學、電學微觀儀器為代表.掃描電子顯微鏡、透射電子顯微鏡(transmission electron microscope,TEM)、原子力顯微鏡等能夠獲取水合物表面微觀結構信息,X 射線衍射(Xray diffraction,XRD)、拉曼譜能夠實現水合物晶體結構參數分析,計算機斷層掃描、核磁共振(nuclear magnetic resonance,NMR) 及成像(magnetic resonance imaging,MRI)則可實現水合物合成及分解過程的結構演化結果,此外也可采用多技術相結合的方式進行綜合分析.

4.1.1 計算機斷層掃描技術

計算機斷層掃描技術被廣泛應用于多學科的微觀結構測試,其基本原理是:基于結構內不同組分對X 或γ 射線的吸收能力不同,通過入射前后射線強度變化的對比分析,獲取結構截面入射方向上的衰減總值,最后利用圖像重建技術得到CT 二維圖像,借助進一步的衰減系數數字矩陣則可得到三維圖像[146].CT 技術的主要優勢在于:(1) 技術成熟、分辨率高,能夠有效開展孔隙尺度的組分及組構分析;(2) 無損觀測,結構完整性高.基于射線掃描,能夠有效保證結構對象的完整性;(3) 結合性高,適用性廣.能夠有效開展基于CT 的聯合實驗,具有廣泛適用性,如CT 三軸剪切實驗[147];(4) 分析結果拓展及延續性強.基于CT 掃描的微觀孔隙結構表征[148]能夠有效用于數值計算等領域研究,如水合物孔隙滲透特性模擬[149].而水合物CT 分析的主要問題為:(1) 水合物與水的區分問題.水合物與水密度相近,因此在CT 圖像中兩者灰度接近,可通過水中添加溶解鹽質[107,150],或使用特殊客體分子[151]實現兩相灰度區分.(2) 低含量水合物識別問題.當水合物含量較低時,其對射線吸收較弱,易被其他結構成分掩蓋.(3) 水合物動態生長、分解過程中的氣水運移及水合物演化識別問題[36].水合物CT 具體應用方面,開展了包括不同水合物孔隙尺度分布[152]、生成及分解過程水合物轉化特征[153-154]、滲透特性孔隙結構關系[149]等相關研究,這些成果為水合物宏觀特性的微觀解釋提供了有力支撐.

4.1.2 掃描電子顯微鏡技術

掃描電子顯微鏡技術通過聚焦電子束對研究對象進行逐點掃描,成像信號可采用二次電子、背散射電子或吸收電子等,并以二次電子為主.利用SEM 可實時獲取研究對象的表面形貌、顆粒關系、孔隙特征等微觀信息.在微觀組構關系研究方面具有十分重要的作用,尤其適用于顆粒及多孔介質研究對象[155].通過超低溫環境實現常壓條件下的水合物穩定存在,形成專用于水合物微觀結構分析的低溫掃描電子顯微鏡(Cryo-SEM),而依據發射電子的不同可進一步分為鎢燈絲型及場發射型掃描電鏡.在實際操作過程中,Cryo-SEM 主要難點在于[36]:(1) 由于超低溫環境,制樣及轉移過程中可能發生水蒸氣表面凝結,需采取干燥空氣等措施避免表面結霜;(2) 噴鍍作業時,注意試樣表面結構完整性,并保證全程超低溫環境;(3) 測試過程避免試樣受電子束的損毀破壞,對電子束操作要求較高.SEM 技術在天然氣水合物賦存特征[156]、多研究角度下的水合物表面結構特征演化[157]、分解動力學分析[158]等方面得到了廣泛應用.

4.1.3 XRD 技術

XRD 技術的原理是基于晶體的周期結構而產生的X 衍射現象,經晶體結構后電磁輻射發生入射方向及強度的變化.被廣泛應用于物質結構的定性和定量分析、微觀結構分類、完整程度分析等,在化學、力學、醫學、生物學、材料學等多學科領域發揮著重要作用.天然氣水合物XRD 技術一般通過液氮實現低溫,能夠有效用于區分水合物晶體類型,是研究水合物晶體結構參數的重要技術手段.水合物XRD 技術實際應用時,要密切關注測試條件度實驗結果的影響,尤其是對XRD 圖譜的影響.重點關注的測試條件包括:步長選擇、掃描速度、累計次數以及測試溫度等.目前XRD 技術在水合物晶體結構類型識別[159]、晶體參數分析[160]、生成分解動力學研究[161]以及“自我保護效應”機理分析[162]等方面取得了重要的成果.除上述介紹外,水合物微觀研究還包括NMR,MRI,拉曼譜等[163-165]微觀測試技術.

微觀測試技術與宏觀力學實驗相結合的宏?微觀研究方式有多種,其中CT 三軸應用最為廣泛.CT 三軸實驗將宏觀三軸剪切實驗與微觀CT 掃描技術相結合,在力學特性分析的基礎上結合剪切過程的微觀結構觀測,實現剪切破壞全過程的宏?微觀分析.

4.2 直接測試技術

與間接測試技術需要結合宏觀力學實驗不同,直接測試技術將微觀測試與微觀力學實驗相統一,能夠直接獲取微觀力學特性及微觀結構演變.相對而言,該方式在直接獲取微觀力學性質方面優勢明顯,而其主要難點在于力學實驗的微觀實現.在微觀化的過程中,不僅伴隨實驗裝置設計及制造精度的提高,實驗操作的難度也將進一步加大.此外一些宏觀實驗原理在微觀層面也將不再適用,這些都極大地限制了直接測試技術的發展,也因此,目前可用于天然氣水合物直接測試的微觀裝置相對較少.

4.2.1 原子力顯微鏡

原子力顯微鏡最早由IBM 在20 世紀末所發明,其基礎是掃描隧道顯微鏡(scanning tunneling microscopy,STM),兩者則被合稱為掃描探針顯微鏡(scanning probe microscopy,SPM)[166].AFM 將微型測力裝置與電子顯微鏡相結合,能有效開展包括絕緣體在內的各種以固體為研究對象的微觀力學實驗,在多學科領域得到了廣泛應用.天然氣水合物研究方面,通過增加溫壓環境,則可開展各種氣體、非氣體水合物的微觀力學特性研究.

AFM 的結構示意圖可見圖12,其主要結構包括微懸臂、激光發射器、激光探測器及載物臺等[167].其中微懸臂為核心部件,其一端固定另一端則有微小的針尖,即探針.探針的運動由控制器通過調節電流實現控制.激光發射器將發射的激光照射在微懸臂的背面,并反射到由光敏二極管所組成的探測器上,這一發射、接收過程即實現了微小應變的光學信號放大.在樣品接觸掃描過程中,探針與樣品表面相互作用,并引起探針位移(彎曲程度)的變化,進而造成反射光束的偏移,光束偏移蘊含的樣品表面信息則被激光探測器所記錄.在這一過程中,控制器會自動調節微懸臂的彎曲程度,保持探針與樣品表面的排斥力恒定,即兩者相對位置始終處于合理范圍內.AFM 的工作模式主要有接觸模式、非接觸模式以及輕敲模式[168](介于接觸模式與非接觸模式之間),不同模式適用于不同的研究對象,并分別具有不同的優缺點,可根據研究需要進行選取.

圖12 AFM 結構示意圖[167]Fig.12 Schematic diagram of AFM[167]

對AFM 進行改造,增加溫壓控制系統,則可形成專用于水合物微觀力學特性研究的水合物AFM(如圖13 所示),當前,部分學者已經開展了相關研究.例如借助于水合物AFM(二氧化硅微球探針),人們發現微球探針的壓入導致了水合物樣品相變的發生,并引起顯著的塑性變形,探針與水合物之間的黏附力主要由水合物表面似液層及分解液形成的液橋所引起[46].當采用AFM 進行水合物的壓痕測量時,應當考慮AFM 尖端引起的熱融化效應[169].Huang等[170]則對THF 水合物結晶過程中的微納米氣泡影響過程進行了研究,重點探討了納米氣泡與水合物間接觸角、納米氣泡周圍水合物晶體分布模式等演化規律.這些研究結論對于進一步分析水合物微觀變形機理、顆粒聚集以及黏附力演化規律等均具有重要的指導意義.

圖13 水合物AFM 結構示意圖[46]Fig.13 Schematic diagram of hydrate AFM[46]

4.2.2 其他測量儀器

除較為常用的AFM 外,一些具有自主研發性質的“微觀機械力學裝置”也被用于天然氣水合物的微觀力學研究.其中,較為廣泛的是基于胡克定律而被研發的一種“微機械測力裝置”[44],其基本的工作原理如圖14 所示.首先在懸臂終端生成水合物顆粒,隨后水合物顆粒發生接觸與分離,在此過程中記錄低彈性常數懸臂終端的顆粒位移,最后根據顆粒位移、懸臂彈性常數,并結合胡克定律則可獲得顆粒之間的黏附力.一般此類實驗需要進行多組實驗,根據實驗結果的平均值進行具體分析[171].此外,一些以上述“微機械測力裝置”為基礎的“改進升級”裝置也被相繼研發[86],在此不再贅述.

圖14 含水合物土的微機械測力裝置原理圖[44]Fig.14 Schematic diagram of devices to measure force in hydrate-bearing soils[44]

基于宏觀三軸實驗原理,日本科研人員研發了一種微三軸實驗裝置[81],可對直徑5 mm,高度10 mm 的試樣開展三軸剪切實驗.同時結合CT 技術,借助于亞毫米及微米尺度圖像,實現剪切面、顆粒變形及旋轉等局部變形的量化分析.Seol 等[84]也提出了一種全新的CT 結合微三軸實驗設備,能夠實現成分替換、流體轉移、顆粒旋轉及破碎的動態觀測.

水合物微觀力學研究相比宏觀性質而言,除了直接的微觀力學關系外,還需要對間接的微觀細節信息(如組構關系、微觀相態分布、動態演化等)進行匯總分析,這些研究內容的有效開展很大程度上依賴于微觀力學裝置(包括微觀測試技術).通過上述總結可知,水合物微觀力學研究領域主要集中在近十年范圍內,而這一時期也是微觀力學裝置快速發展的階段.在宏微觀相結合的趨勢下,實驗裝置匱乏與微觀力學發展訴求之間的不協調應當引起廣大學者們的重視,并作為下一階段重點解決的問題之一.

5 含水合物土微觀力學研究挑戰

圍繞含水合物土微觀力學特性,廣大學者開展了從原子(分子)到晶體,從晶體到顆粒(界面),從顆粒到整體(含水合物土)的系統性研究,在取得豐碩研究成果的同時,也面臨著諸多挑戰.

(1)含水合物土微觀結構辨識及定量表征

含水合物土的內部結構性顯著,其宏觀力學行為響應很大程度上受控于微觀結構演化,微觀結構的有效測量、準確辨識及定量表征是微觀力學特性研究的基礎.然而,目前含水合物土微觀力學特性研究更多是基于孔隙尺度的內部結構測量、辨識與表征,但是更小尺度的內部結構數據匱乏,導致純水合物及其與土顆粒界面的力學行為仍不甚清楚.

(2)非均質含水合物土微觀力學特性研究

實際自然環境中,水合物生成受多種因素的影響,含水合物土呈現出“非均質性”,決定了其力學特性也是各向異性的.由于制樣方法或實驗條件的不足,現階段針對含水合物土的微觀力學特性研究,幾乎都是以均質水合物為研究對象,如基于均勻制樣的力學實驗及均勻成核的MD 研究等,迫切需要深入開展非均質含水合物土力學特性研究.

(3)含水合物土多場耦合微觀力學研究

實際工程環境中,水合物開采是一個多場耦合的過程,儲層變形場與其他物理場之間的耦合作用規律是實現工程地質風險可控的關鍵.然而,考慮多場耦合作用的含水合物土微觀力學實驗研究涉及較少,目前僅通過數值模擬的手段實現了力學多場耦合分析,相應的微觀力學實驗仍然是一個極具挑戰性的課題.

(4)微觀本構模型構建及拓展研究

本構模型方面,目前已經建立了部分考慮水合物微觀孔隙結構以及反映微觀力學特性的本構模型,而其適用范圍以及參數的確定則是本構研究的主要難點,當前仍缺乏具有廣泛適用性同時具有有限模型參數的微觀本構模型.另一方面,圍繞本構模型的系統性研究不足,基于本構模型的二次開發及其數值模擬研究較為鮮見,這無疑限制了本構模型的進一步應用與推廣.

6 結論與建議

天然氣水合物的開采過程是多相多物理場的動態演化過程,在相變、滲流、傳熱等多因素作用下,微觀孔隙結構持續變化,最終引起宏觀力學特性的演變.因此,天然氣水合物宏觀力學特性的內在本質是其微觀組構關系,微觀力學特性研究的重要性不言而喻.整體而言,在宏觀力學研究的基礎上,水合物的微觀力學研究已經進入快車道,同時在宏微觀相結合的趨勢下,微觀、介觀到宏觀的力學特性統一問題已經成為廣大學者研究的焦點之一.而多學科相交叉的綜合研究也將進一步推動天然氣水合物微觀力學研究的快速發展.在本文系統總結天然氣水合物微觀力學研究現狀的基礎上,并針對目前研究存在的不足,得到如下幾點認識或建議:

(1)微觀力學特性測試裝置及技術.含水合物土微觀力學研究對實驗裝置及其技術的依賴度較高,目前可用于直接測量微觀力學特性的實驗裝置相對匱乏.此外在多學科交叉背景下,相關學科的最新微觀測試技術尤其是光、電、聲學測量儀器,應盡快應用到水合物微觀力學研究領域.因此,在對現有實驗裝置進行升級與改造的基礎上,應當大力開展水合物微觀力學實驗裝置的研發工作,這也是含水合物土微觀力學發展的必由之路.

(2)多種研究方式相結合.含水合物土微觀力學研究是系統性體系研究,然而受各種因素及實際操作的限制,不同研究方法之間尚缺乏有效聯系的紐帶,不同研究方法所得結論相對獨立并缺乏對比性.而不同研究方法之間的成果互通無疑是一種有效的研究捷徑,同時彼此之間的結果印證也將使得研究結論更具系統性.例如非均勻制樣及特殊儲層制樣是實驗研究的一大難題,而根據PFC 研究,通過合理的模型建立及參數選取,能夠有效解決制樣難題,在定量表征方面也具有一定優勢;基于理論分析的本構模型建立可對實驗結果進行修正與完善,其二次開發又可實現數值模擬分析.

(3)微觀力學特性本構研究.在構建反映微觀力學特性本構模型方面,應當避免過多的模型參數.首先對建模框架進行優選,選取相對成熟并且具有廣泛適用性的經典本構模型作為建模框架,通過引入微觀關聯參數構建起反應微觀特性的本構模型.而模型參數應盡可能具有實際物理意義,盡量避免“硬性”參數的引入,降低參數確定難度.此外還應開展基于微觀本構模型的優選方法研究,尤其是針對不同儲層特性、賦存模式以及多因素影響,這需要進一步開展大量的微觀力學實驗研究.

(4)典型儲層及賦存模式研究.與現階段所取得的微觀力學研究成果相比,典型水合物儲層研究略顯不足.尤其是全球水合物資源的90%以上賦存于泥質粉砂型儲層(裂隙型水合物)中,而微觀力學研究尚未開展針對泥質粉砂水合物的有效研究,因此應逐步開展泥質粉砂等典型水合物儲層的微觀力學研究.特別是針對開采過程中的水合物賦存模式動態演化微觀力學研究也應當盡快開展.針對上述研究可重點關注DEM 及微觀測試技術的應用.

猜你喜歡
實驗模型研究
一半模型
記一次有趣的實驗
FMS與YBT相關性的實證研究
遼代千人邑研究述論
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
做個怪怪長實驗
EMA伺服控制系統研究
3D打印中的模型分割與打包
主站蜘蛛池模板: 午夜老司机永久免费看片| 狠狠综合久久久久综| 最新国产你懂的在线网址| 午夜影院a级片| 欧美色99| 青青操国产视频| 日本手机在线视频| 国产精品19p| 亚洲国产精品国自产拍A| AV不卡在线永久免费观看| 亚洲综合天堂网| a级毛片免费网站| 亚洲色中色| 一本久道久综合久久鬼色| 亚洲色婷婷一区二区| 欧美日韩精品综合在线一区| 亚洲视频免费在线看| 丝袜亚洲综合| 国产人在线成免费视频| 中文字幕调教一区二区视频| 国产免费怡红院视频| 亚洲国产欧美国产综合久久 | 久久先锋资源| 欧美激情伊人| 国产欧美成人不卡视频| a在线观看免费| 久久精品国产亚洲麻豆| 欧美中文字幕无线码视频| 99久久性生片| 97在线视频免费观看| 亚洲天堂视频在线播放| 成人无码一区二区三区视频在线观看| www.亚洲天堂| 幺女国产一级毛片| 国产精品亚洲专区一区| 国产精品欧美在线观看| 波多野结衣无码AV在线| 国产女人爽到高潮的免费视频 | 一级毛片视频免费| 国产不卡网| 国产精品yjizz视频网一二区| 国产成人欧美| 欧美午夜一区| 一区二区午夜| 中文字幕第1页在线播| 国内精品一区二区在线观看| 欧美日韩国产成人高清视频| 亚洲精品在线91| 色网站在线视频| www亚洲天堂| 国产美女一级毛片| 亚洲一区网站| 国产精品尤物在线| 无码精品福利一区二区三区| 再看日本中文字幕在线观看| 久久精品嫩草研究院| 国产欧美日韩综合在线第一| 亚洲人成网线在线播放va| 国产成人啪视频一区二区三区| 在线国产资源| AV熟女乱| 亚洲精品国产成人7777| 亚洲欧美日韩天堂| 久久精品丝袜高跟鞋| 一级一级特黄女人精品毛片| 久久国产黑丝袜视频| 国产在线麻豆波多野结衣| 四虎永久在线精品国产免费| 老色鬼欧美精品| 特级aaaaaaaaa毛片免费视频 | 日韩第一页在线| 国产精品丝袜在线| 麻豆精品在线视频| 亚洲欧洲日韩综合| 中国黄色一级视频| 国产精品自在自线免费观看| 精品久久久久久中文字幕女| 亚洲日韩AV无码一区二区三区人 | 国产精品亚洲а∨天堂免下载| 久久精品丝袜| 亚洲第一视频免费在线| 亚洲国产91人成在线|