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

    基于離散元的擠壓構(gòu)造定量分析與模擬

    2022-08-31 02:51:14李長(zhǎng)圣尹宏偉徐雯嶠吳珍云管樹(shù)巍
    大地構(gòu)造與成礦學(xué) 2022年4期
    關(guān)鍵詞:變形模型

    李長(zhǎng)圣, 尹宏偉, 徐雯嶠, 吳珍云, 管樹(shù)巍, 賈 東, 任 榮

    基于離散元的擠壓構(gòu)造定量分析與模擬

    李長(zhǎng)圣1, 2, 3, 4, 尹宏偉3*, 徐雯嶠3, 吳珍云1, 2, 管樹(shù)巍4, 賈 東3, 任 榮4

    (1. 核資源與環(huán)境國(guó)家重點(diǎn)實(shí)驗(yàn)室, 江西 南昌 330013; 2. 東華理工大學(xué)地球科學(xué)學(xué)院, 江西 南昌 330013; 3. 南京大學(xué) 地球科學(xué)與工程學(xué)院, 江蘇 南京 210023; 4. 中國(guó)石油勘探開(kāi)發(fā)研究院, 北京 100083)

    隨著離散元理論和計(jì)算機(jī)技術(shù)的發(fā)展, 離散元方法已經(jīng)廣泛應(yīng)用到不同尺度的構(gòu)造模擬中。相較于傳統(tǒng)的沙箱模擬實(shí)驗(yàn), 離散元可以更精確地控制實(shí)驗(yàn)的邊界條件, 定量的分析構(gòu)造變形過(guò)程, 有助于從細(xì)觀尺度深入認(rèn)識(shí)地層力學(xué)性質(zhì)對(duì)構(gòu)造變形過(guò)程的影響。本文系統(tǒng)闡述了基于離散元的構(gòu)造變形定量分析方法, 結(jié)合一個(gè)典型的擠壓構(gòu)造離散元數(shù)值模擬實(shí)驗(yàn), 模擬了水平擠壓環(huán)境下構(gòu)造的形成過(guò)程, 并對(duì)變形過(guò)程中的應(yīng)力、應(yīng)變分布變化與裂縫生成規(guī)律進(jìn)行分析, 取得以下認(rèn)識(shí): (1) 該模擬實(shí)驗(yàn)的構(gòu)造以前展式的逆沖疊瓦斷層為主, 斷層從擠壓端到遠(yuǎn)離擠壓端依次活動(dòng); (2) 裂縫與斷層形成有密切關(guān)系, 局部區(qū)域內(nèi)聚集的大量裂縫是產(chǎn)生斷層的誘因; (3) 斷層形成初期,斷層內(nèi)物質(zhì)運(yùn)動(dòng)距離較小, 產(chǎn)生裂縫增量最多; 斷層活動(dòng)后期, 裂縫增量減少; (4) 體積應(yīng)變可以表征裂縫類型(拉張或壓縮), 變形應(yīng)變可以區(qū)分正向和反向逆沖斷層; (5) 平均應(yīng)力大小與地形起伏呈正相關(guān)關(guān)系, 最大剪切應(yīng)力持續(xù)在將要形成的新斷層處累積, 直至該新斷層形成, 最大剪切應(yīng)力繼而消散, 繼續(xù)往前傳播, 在下一個(gè)將要形成的新斷層處累積。研究結(jié)果表明離散元方法在構(gòu)造變形、應(yīng)力應(yīng)變與裂縫預(yù)測(cè)定量研究中具有巨大潛力。

    離散元; 構(gòu)造變形; 擠壓構(gòu)造; 應(yīng)力應(yīng)變; 定量分析

    0 引 言

    離散元方法(discrete element method, 簡(jiǎn)稱DEM)是Cundall and Strack (1979)提出的用于研究巖土體變形的一種方法, 其基本思想是把自然界的巖土體看成是離散的單元幾何體。離散元方法可以有效地模擬沉積地層中出現(xiàn)的斷層及斷層相關(guān)褶皺等脆性變形(Dean et al., 2013; 李長(zhǎng)圣, 2019), 廣泛應(yīng)用于解決構(gòu)造相關(guān)的地質(zhì)問(wèn)題。例如, 分析鹽刺穿過(guò)程中沉積蓋層的破裂機(jī)制(張潔等, 2008; Yin et al., 2009), 研究地層內(nèi)聚力對(duì)構(gòu)造形態(tài)及應(yīng)力應(yīng)變分布的影響(Morgan, 2015), 揭示水平擠壓環(huán)境下滑脫層對(duì)構(gòu)造變形過(guò)程中的應(yīng)變分布變化與裂縫生成規(guī)律(蔡申陽(yáng)等, 2016), 解析裂陷盆地演化過(guò)程中分層伸展疊加變形的動(dòng)力學(xué)演化過(guò)程(周易等, 2019), 分析區(qū)域應(yīng)力場(chǎng)作用下反轉(zhuǎn)構(gòu)造特性及正斷層反轉(zhuǎn)控震機(jī)制(吳珍云等, 2019), 揭示純走滑拉分盆地發(fā)育過(guò)程中的斷裂擴(kuò)展和連接過(guò)程(劉源和Konietzky, 2019), 探討庫(kù)車前陸沖斷帶西部鹽構(gòu)造形成的控制因素及其形成機(jī)理(李維波等, 2017; 李江海等, 2020; 徐雯嶠等, 2020)。

    Schreurs et al. (2006, 2016)通過(guò)對(duì)比全球多個(gè)物理模擬實(shí)驗(yàn)室的擠壓構(gòu)造實(shí)驗(yàn), 發(fā)現(xiàn)在相同的實(shí)驗(yàn)條件和方法下, 不同的實(shí)驗(yàn)?zāi)M結(jié)果無(wú)法完全一致。而DEM在相同的初始條件(顆粒位置和半徑相同)和邊界條件下, 實(shí)驗(yàn)完全可重復(fù)。其中, DEM顆粒材料屬性通過(guò)細(xì)觀參數(shù)標(biāo)定, 可選擇的材料較多。并且, 所有的變量(如位移、應(yīng)力、應(yīng)變、速度等信息)都可以實(shí)時(shí)監(jiān)測(cè)。在物理模擬中, 變形需要通過(guò)圖像分析(如粒子圖像測(cè)速)(沈禮等, 2012)、激光掃描(微型激光測(cè)高)(Liu et al., 2013)、利用計(jì)算機(jī)X射線斷層掃描技術(shù)進(jìn)行體掃描(李長(zhǎng)圣, 2014; 李長(zhǎng)圣等, 2014; Li et al., 2016)等方法獲取。如果沒(méi)有這些設(shè)備, 則只能通過(guò)對(duì)模型側(cè)面拍照來(lái)觀察模型側(cè)面形變, 或切割模型觀察其內(nèi)部形變。而DEM則可以給出每個(gè)變形階段的所有信息, 用這些信息可以計(jì)算出系統(tǒng)的應(yīng)力應(yīng)變場(chǎng)。同時(shí), 同一初始模型, 模擬結(jié)果一致, 利于單因素變量(如擠壓速率、古隆起、先存斷層、同構(gòu)造沉積與剝蝕等)研究(李長(zhǎng)圣, 2019; 辛文等, 2020; 徐雯嶠等, 2020; Li et al., 2021; Xu et al., 2021; 鄒瑋等, 2021)。

    離散元數(shù)值模擬方法在構(gòu)造形態(tài)、斷裂預(yù)測(cè)及應(yīng)力應(yīng)變的定量分析方面具有明顯的優(yōu)勢(shì)。眾多學(xué)者已經(jīng)將離散元引入構(gòu)造變形的模擬中, 但分析的重點(diǎn)多集中于構(gòu)造幾何形態(tài)的定性解析。本文在闡述離散元原理、模擬過(guò)程的基礎(chǔ)上, 結(jié)合一個(gè)典型的擠壓構(gòu)造模擬實(shí)例, 詳細(xì)剖析了基于離散元的擠壓構(gòu)造實(shí)驗(yàn)過(guò)程及定量分析方法。

    1 基本原理

    DEM基本思想是將顆粒材料內(nèi)部細(xì)觀尺度的單個(gè)離散顆粒視為一個(gè)離散單元, 通過(guò)一系列離散的單元來(lái)模擬顆粒材料的力學(xué)行為。離散元的求解實(shí)際上是迭代計(jì)算顆粒位移和受力, 可以概括為兩部分: 第一步, 已知的顆粒所受合力和合力矩, 由牛頓第二定律更新每個(gè)顆粒的位置; 第二步, 找到相互接觸的顆粒, 應(yīng)用接觸力學(xué)模型(即力?位移法則)。本文采用Hertz-Mindlin接觸力學(xué)模型, 其在力的計(jì)算方面精確且高效, 同時(shí)引入了粘結(jié)(Bond)接觸力學(xué)模型(簡(jiǎn)稱粘結(jié)模型), 以模型脆性巖石的變形行為, 模擬效果較好。

    相互粘結(jié)的兩個(gè)顆粒, 其接觸力的計(jì)算可以分為兩種狀態(tài), 即壓縮與拉伸。當(dāng)兩個(gè)顆粒處于壓縮狀態(tài)下, 接觸力計(jì)算采用Hertz-Mindlin接觸力學(xué)模型; 當(dāng)兩個(gè)顆粒處于拉伸狀態(tài)下, 接觸力計(jì)算采用粘結(jié)模型。

    一般地, 規(guī)定兩個(gè)顆粒間壓力為正, 拉力為負(fù)。

    式中:U. 顆粒間的法向疊合量;K. 顆粒間的法向接觸剛度。

    式中:U. 顆粒間的切向疊合量;K. 顆粒間的切向接觸剛度。

    (1) 壓縮狀態(tài)

    當(dāng)兩個(gè)顆粒(顆粒O和顆粒A, 圖1)處于壓縮狀態(tài)下時(shí),K由下式計(jì)算,

    其中,

    式中:O和A分別為顆粒O和A的剪切模量,和分別為顆粒O和A的泊松比。

    當(dāng)顆粒處于壓縮狀態(tài)下時(shí),K由下式計(jì)算,

    (2) 拉伸狀態(tài)

    相互粘結(jié)的兩個(gè)顆粒處于拉伸狀態(tài)下時(shí), 法向接觸剛度K和切向接觸剛度K采用以下公式計(jì)算,

    圖1 接觸面定義

    式中:E. 粘結(jié)的楊氏模量;G. 粘結(jié)的聚合模量;. 兩個(gè)接觸顆粒的接觸面的面積, 為一個(gè)圓, 該圓半徑簡(jiǎn)化為較小顆粒的半徑o, 如圖1所示,o2;o和A分別為顆粒O和A的半徑。

    當(dāng)顆粒處于粘結(jié)狀態(tài)時(shí), 粘結(jié)的抗拉力為,

    粘結(jié)的抗剪力為,

    式中:T.粘結(jié)的抗拉強(qiáng)度;C.粘結(jié)的剪切強(qiáng)度。相互粘結(jié)的顆粒, 自始至終滿足FFmax并且|F|≤Fmax時(shí), 顆粒處于粘結(jié)狀態(tài)。當(dāng)上述一個(gè)條件不滿足時(shí), 顆粒間的粘結(jié)斷開(kāi), 進(jìn)入非粘結(jié)狀態(tài), 之后不再重新產(chǎn)生粘結(jié)。

    當(dāng)兩個(gè)顆粒處于非粘結(jié)狀態(tài)(斷裂)時(shí),Fmax= 0.0, 顆粒間不產(chǎn)生拉力;Fmax=μF, 切向力需滿足|F|≤Fmax。

    本文的數(shù)值模擬使用李長(zhǎng)圣開(kāi)發(fā)的離散元軟件(李長(zhǎng)圣等, 2017; Li et al., 2018, 2021; 李長(zhǎng)圣, 2019), 選用Hertz-Mindlin接觸力學(xué)模型和粘結(jié)(Bond)模型, 計(jì)算顆粒間的接觸力(Morgan, 2015; 李長(zhǎng)圣, 2019), 采用蛙跳法更新顆粒位移(Itasca Consulting Group, 2008; 李長(zhǎng)圣, 2019)。

    2 應(yīng)力應(yīng)變表征

    應(yīng)力應(yīng)變作為描述材料受力和變形的基本力學(xué)量, 在連續(xù)介質(zhì)力學(xué)中被廣泛應(yīng)用, 但由于顆粒材料的離散性, 基于連續(xù)介質(zhì)定義的應(yīng)力應(yīng)變無(wú)法直接描述離散顆粒集合體的本構(gòu)特征。本文采用一種在構(gòu)造模擬中應(yīng)用效果較為理想的等效應(yīng)力和等效應(yīng)變定義方法(劉其鵬, 2010)。如沒(méi)有特殊說(shuō)明, 本文提到的應(yīng)力應(yīng)變即為等效應(yīng)力和等效應(yīng)變。

    2.1 應(yīng) 變

    連續(xù)介質(zhì)力學(xué)中, 應(yīng)變是局部變形區(qū)域的位移梯度(O’Sullivan, 2011)。DEM中, 可以通過(guò)追蹤顆粒每一時(shí)步的位置, 得到顆粒的累積位移, 并分成水平位移和垂直位移兩部分。采用臨近點(diǎn)插值算法(MathWorks, 2015), 將位移投射到規(guī)則的笛卡爾網(wǎng)格中。構(gòu)造研究中, 地層變形一般較大, 需采用有限應(yīng)變理論進(jìn)行計(jì)算分析(Oertel, 1996; Means, 2012)。

    集合體隨時(shí)間發(fā)生變化, 可以分解為體積變化導(dǎo)致的形變(體積應(yīng)變)和剪切變形引起的形變(剪切應(yīng)變), 本文采用Morgan (2015)給出的方法計(jì)算應(yīng)變。把每個(gè)顆粒位移分成水平和垂直兩部分, 其二維的位移梯度張量為D,D=?U/?U(指標(biāo)符號(hào),= 1, 2), 其中U為位移矢量網(wǎng)格。有限應(yīng)變張量的不變量用來(lái)表征變形場(chǎng)在時(shí)間和空間上的變化。有限應(yīng)變的第一不變量I是體積應(yīng)變(Malvern, 1969; Jaeger et al., 2009), 由位移梯度張量得出。F=?x/?X,xX分別為變形后和變形前的顆粒坐標(biāo)。的行列式形式為,

    它表征了集合體體積的變化率。

    體積應(yīng)變I計(jì)算方法如下:

    變形應(yīng)變用有限應(yīng)變張量的第二不變量II′表征, 其中′=–I/2。II′由的分量按下式計(jì)算,

    本文采用體積應(yīng)變I(表征體積變化)和有限應(yīng)變張量的第二不變量II′(表征剪切形變程度)描述計(jì)算結(jié)果。

    2.2 應(yīng) 力

    考慮一個(gè)由圓形顆粒組成的體積為(包括孔隙)的顆粒集合即表征元, 將顆粒集合視為連續(xù)體, 粒材料的應(yīng)力定義為集合內(nèi)的平均應(yīng)力, 假設(shè)個(gè)接觸作用在個(gè)顆粒上, 則(Morgan, 2015; 李長(zhǎng)圣, 2019)

    式中:σ. 最大主應(yīng)力;σ. 最小主應(yīng)力; 最大剪切應(yīng)力為max=Δ/2。

    式中:r. 顆粒的半徑大小;. 集合體的顆粒數(shù)量。

    3 實(shí)驗(yàn)過(guò)程

    3.1 參數(shù)標(biāo)定

    與有限元、有限差分等連續(xù)介質(zhì)力學(xué)方法不同, 離散元中采用顆粒的細(xì)觀參數(shù)(如顆粒間的剛度系數(shù)、摩擦系數(shù)等)表征顆粒系統(tǒng)的宏觀力學(xué)性質(zhì)。為了模擬自然界地層的構(gòu)造變形行為, 通過(guò)雙軸實(shí)驗(yàn)得到離散元顆粒的細(xì)觀參數(shù)與離散元顆粒體所表現(xiàn)出的宏觀響應(yīng)間的關(guān)系, 即細(xì)觀參數(shù)(如顆粒間的摩擦系數(shù)、剪切模量)與宏觀參數(shù)(如粘聚力、內(nèi)摩擦角等)間的映射關(guān)系。通過(guò)5組雙軸試驗(yàn), 得到如表2所示的細(xì)觀參數(shù)和宏觀參數(shù)對(duì)應(yīng)關(guān)系, 詳細(xì)的測(cè)試過(guò)程和結(jié)果見(jiàn)李長(zhǎng)圣(2019)。雙軸實(shí)驗(yàn)的應(yīng)力應(yīng)變及剪切強(qiáng)度包絡(luò)線見(jiàn)圖2, 顆粒的細(xì)觀參數(shù)見(jiàn)表1, 接觸的粘結(jié)(Bond)參數(shù)見(jiàn)表2。

    3.2 擠壓實(shí)驗(yàn)

    (1) 在長(zhǎng)60 km, 高15 km的矩形盒子內(nèi), 由92個(gè)直徑160 m的顆粒分別生成左墻和右墻,由749個(gè)直徑80 m顆粒生成底墻。之后, 在該矩形盒子內(nèi), 按照1∶1的比例隨機(jī)生成直徑為120 m和160 m的兩種大小顆粒。新生成顆粒與之前已經(jīng)生成的顆粒進(jìn)行疊合判斷, 如果新生成的顆粒與已生成的顆粒疊合, 則刪除該顆粒, 重新生成一個(gè)顆粒。當(dāng)顆粒生成疊合判斷的嘗試次數(shù)大于2000次(Itasca Consulting Group, 2008; 李長(zhǎng)圣, 2019), 則停止生成顆粒, 共生成25027個(gè)顆粒。

    (2) 設(shè)置顆粒間的細(xì)觀參數(shù), 具體值見(jiàn)表1, 所有顆粒的摩擦系數(shù)設(shè)置為0.0。

    (3) 顆粒生成之后, 將在重力作用下做自由落體運(yùn)動(dòng), 同時(shí)在阻尼力的作用下, 約5萬(wàn)步后, 沉積穩(wěn)定。之后, 刪除縱坐標(biāo)5 km以上的顆粒(不刪除墻體), 體系上覆地層重量減小, 體系有微小的回彈膨脹, 繼續(xù)計(jì)算1萬(wàn)步, 保證體系穩(wěn)定。最終, 得到長(zhǎng)60 km, 高5 km的初始模型(圖3), 共18606個(gè)顆粒。將左右邊墻和巖層顆粒摩擦系數(shù)設(shè)為0.3, 巖層顆粒細(xì)觀參數(shù)見(jiàn)表1、2。同時(shí), 設(shè)置顆粒顏色, 以便于區(qū)分巖層。

    (4) 底墻和右墻固定, 給定左邊墻2 m/s水平向右的恒定速度, 當(dāng)左邊墻向右移動(dòng)20 km時(shí), 結(jié)束計(jì)算。

    4 結(jié)果分析

    擠壓實(shí)驗(yàn)的構(gòu)造演化過(guò)程見(jiàn)圖4。隨著擠壓的進(jìn)行, 斷層F1、F2、F3和F4依次在模型收縮1 km、4 km、9 km和16 km后開(kāi)始強(qiáng)烈活動(dòng), 產(chǎn)生明顯的斷距, 最終形成典型的疊瓦狀逆沖斷層帶, 該斷層帶整體呈現(xiàn)出由擠壓端向模型內(nèi)部高度逐漸減薄的楔體構(gòu)造形態(tài), 且楔體構(gòu)造外無(wú)明顯的構(gòu)造變形(圖6e)。受擠壓端邊界效應(yīng)的影響, 在靠近擠壓端產(chǎn)生一條反沖斷層, 與物理模擬結(jié)果一致(Sun et al., 2016; 李長(zhǎng)圣, 2019; Cui et al., 2020; Li, 2021)。

    圖2 雙軸實(shí)驗(yàn)結(jié)果

    表1 顆粒的細(xì)觀參數(shù)

    注:按照1∶1的比例隨機(jī)生成直徑為120 m和160 m的兩種大小顆粒。

    表2 顆粒間的粘結(jié)參數(shù)

    注:壓縮速度取2.0 m/s, 實(shí)驗(yàn)過(guò)程是準(zhǔn)靜態(tài)過(guò)程, 詳細(xì)調(diào)試過(guò)程及結(jié)果見(jiàn)李長(zhǎng)圣(2019)。

    圖3 初始模型

    圖4 模型收縮量為1 km (a)、4 km (b)、9 km (c)、16 km (d)、20 km (e)時(shí),楔體的構(gòu)造解釋

    4.1 楔體構(gòu)造定量分析

    楔體的組成要素包括坡頂、坡趾、坡角和坡面(圖4c)。離散元中, 研究對(duì)象為離散顆粒, 可采用基于網(wǎng)格的搜尋地表法識(shí)別出坡面, 采用最遠(yuǎn)距離法搜尋坡趾。與前人研究(Buiter et al., 2006, 2016; Schreurs et al., 2006; Schreurs, 2016; Sun et al., 2016; Wang et al., 2021)相比, 該方法給出了坡趾嚴(yán)格定義, 適合編程實(shí)現(xiàn)自動(dòng)測(cè)量, 排除了人為因素導(dǎo)致的坡角測(cè)量誤差, 詳細(xì)描述見(jiàn)李長(zhǎng)圣(2019)和Li et al. (2021)。

    由楔體的高度(圖5a)和寬度(圖5b)變化可知, 楔體高度前期呈線性增加, 后期趨于穩(wěn)定; 楔體寬度則呈現(xiàn)階梯式的增加, 在斷層形成(收縮量1 km、4 km、9 km和16 km)時(shí), 楔體寬度陡增, 物理模擬和數(shù)值模擬均表明楔體寬度存在這種周期性變化的現(xiàn)象(Bigi et al., 2010; Buiter et al., 2016; Schreurs et al., 2016; Sun et al., 2016; Li et al., 2021; Wang et al., 2021), 并且與斷層形成存在密切關(guān)系。與楔體高度和寬度不同, 楔體坡角較快進(jìn)入一個(gè)穩(wěn)定值(圖5c), 總體呈現(xiàn)波動(dòng)狀態(tài), 這種楔體坡角波動(dòng)可以總結(jié)為兩個(gè)狀態(tài): 首先隨著擠壓進(jìn)行楔體逐漸增厚, 坡角顯著增大, 如模型從1 km擠壓到 4 km; 然后, 新的斷層F2形成, 坡角隨即減小, 即模型從4 km收縮到 6 km, 接著楔體繼續(xù)增厚, 坡角顯著增大, 進(jìn)入下一個(gè)演化循環(huán)。新斷層開(kāi)始活動(dòng), 楔體高度持續(xù)增大, 楔體寬度相對(duì)穩(wěn)定, 坡角減小。應(yīng)力累計(jì)、斷裂持續(xù)增加, 形成斷層, 使得坡角發(fā)生波動(dòng), 這種波動(dòng)值大小可能與地層強(qiáng)度有關(guān)。總體上, 坡角趨于一個(gè)穩(wěn)定值, 符合臨界楔理論(Davis et al., 1983)。臨界楔理論已經(jīng)被很多數(shù)值模擬和物理模擬結(jié)果所驗(yàn)證(Buiter, 2012), 成功解釋了如臺(tái)灣造山帶等活動(dòng)邊緣中的許多地質(zhì)現(xiàn)象(Davis et al., 1983; Dahlen, 1984; Dahlen, 1990; Suppe, 2007)。如果沒(méi)有新的物質(zhì)介入, 楔體將沿著滑脫層滑動(dòng)而不產(chǎn)生內(nèi)部形變。而處于次臨界和超臨界狀態(tài)的楔體是不穩(wěn)定的, 其會(huì)調(diào)整內(nèi)部形變來(lái)達(dá)到穩(wěn)定狀態(tài)。次臨界楔體通過(guò)產(chǎn)生新的逆沖斷層或者使老的斷層活化, 而使坡角增加; 超臨界楔體通過(guò)前陸形變, 使坡角減小(孫闖, 2017)。所以, 楔體是以局部滑動(dòng)和產(chǎn)生新的逆沖或使老斷層活化抬升的形式周期性演化的(Morgan, 2015; Sun et al., 2016; 李長(zhǎng)圣, 2019; Li et al., 2021; Wang et al., 2021)。事實(shí)上, 楔體幾何形態(tài)是隨時(shí)間一直變化的, 臨界角只是一個(gè)近似值(Gutscher et al., 1996, 1998)。本次實(shí)驗(yàn)中, 模型的臨界楔角約18° (圖5c), 與物理模擬中石英砂的擠壓實(shí)驗(yàn)相近(李長(zhǎng)圣, 2019; Li et al., 2021)。

    圖5 楔體高度(a)、寬度(b)和坡角(c)隨收縮量的變化

    4.2 斷裂解譯

    擠壓過(guò)程中, 楔體的構(gòu)造以前展式的逆沖疊瓦斷層為主(圖4), 靠近擠壓端的斷層到遠(yuǎn)離擠壓端的斷層依次活動(dòng)。由于顆粒間設(shè)置了粘結(jié)(Bond)力, 相互粘結(jié)的顆粒間具有粘結(jié)力鏈, 使得楔體具有抗剪強(qiáng)度和抗拉強(qiáng)度, 斷層形成時(shí), 粘結(jié)力鏈會(huì)發(fā)生剪切或者拉伸斷裂。

    收縮量為1 km、4 km、9 km、16 km、20 km時(shí), 粘結(jié)力鏈分布圖見(jiàn)圖6。模型初始化階段, 模型中相互接觸的顆粒生成粘結(jié)力鏈, 粘結(jié)力鏈數(shù)最大。隨著模型收縮, 粘結(jié)力鏈達(dá)到設(shè)置的抗剪或者抗拉強(qiáng)度, 將有粘結(jié)力鏈斷開(kāi)。顆粒斷開(kāi)后, 可能再次接觸, 但是它們之間不再產(chǎn)生粘結(jié), 這時(shí)它們?yōu)闊o(wú)粘結(jié)接觸。隨著擠壓進(jìn)行, 粘結(jié)力鏈逐漸斷開(kāi), 依次形成斷層F1、F2、F3和F4, 這些斷層處, 地層裂縫發(fā)育密集, 斷層間地層裂縫相對(duì)稀疏。模型收縮量達(dá)到20 km時(shí), F4斷層前端形成一對(duì)“V”字型正反逆沖斷層, 表明斷層并不是從底墻某一破裂點(diǎn)發(fā)育而來(lái), 而是從淺地表發(fā)育密集小斷裂, 并逐漸先下延伸貫穿為斷層, 與Wang et al. (2021)研究結(jié)果不同。

    圖6 收縮量為1 km (a)、4 km (b)、9 km (c)、16 km (d)、20 km (e)時(shí), 粘結(jié)力鏈分布(藍(lán)色為粘結(jié)力鏈, 黃色為無(wú)粘結(jié)的接觸)

    每擠壓1 km為一個(gè)粘結(jié)力鏈數(shù)統(tǒng)計(jì)單位, 每擠壓1 km, 新斷開(kāi)的粘結(jié)力鏈數(shù)見(jiàn)圖7。例如, 模型從0 km收縮到1 km, 新斷開(kāi)的粘結(jié)力鏈數(shù)為1611個(gè), 從3 km到4 km新斷開(kāi)的粘結(jié)力鏈數(shù)為1826個(gè)。模型從0 km收縮到1 km, 斷層F1開(kāi)始活動(dòng); 從3 km收縮到4 km, 斷層F2開(kāi)始活動(dòng); 從8 km收縮到9 km, 斷層F3開(kāi)始活動(dòng); 從15 km收縮到16 km, 斷層F4開(kāi)始活動(dòng)。新斷層的形成具有瞬時(shí)性, 以圖7中F3斷層為例, 新斷層(F3)開(kāi)始活動(dòng)時(shí), 新斷開(kāi)的粘結(jié)力鏈最大, 之后, 新斷開(kāi)的粘結(jié)力鏈數(shù)逐漸減少。也就是說(shuō), 隨著模型逐漸收縮, F3斷層處應(yīng)力累積導(dǎo)致F3形成, 形成密集地層裂縫, 釋放大量能量。之后, 該斷層將逐漸緩慢釋放能量, 直到更新的斷層F4形成, 此時(shí)老斷層F3停止活動(dòng)。

    模型收縮過(guò)程中, 斷層的滑移曲線見(jiàn)圖8。隨著模型收縮, 斷層F1、F2、F3和F4依次產(chǎn)生, 新斷層產(chǎn)生, 老斷層活動(dòng)逐漸停止。模型在收縮到3 km時(shí), 斷層F2可能已經(jīng)開(kāi)始活動(dòng), 模型收縮到4 km時(shí), F2斷層在變形、力鏈和應(yīng)變圖上可見(jiàn); F2斷層形成時(shí), F1斷層仍在活動(dòng), 相鄰斷層活動(dòng)時(shí)間有疊合(圖8灰色部分), 新斷層產(chǎn)生, 老斷層逐漸趨于穩(wěn)定。F3和F2、F4和F3活動(dòng)時(shí)間也有同樣的現(xiàn)象。

    選取斷層F1、F2、F3和F4處的四個(gè)監(jiān)測(cè)點(diǎn)PF1、PF2、PF3和PF4, 繪制了它們的運(yùn)動(dòng)路徑(圖9a)。隨著模型收縮, 點(diǎn)PF1累計(jì)位移量最大, PF2、PF3和PF4依次減小(圖9b)。為了進(jìn)一步分析斷層內(nèi)粘結(jié)力鏈演化過(guò)程, 統(tǒng)計(jì)了各個(gè)監(jiān)測(cè)點(diǎn)周圍2 km范圍內(nèi)的粘結(jié)力鏈斷開(kāi)情況(圖9c)。各個(gè)斷層活動(dòng)時(shí), 其內(nèi)部地層發(fā)育大量斷裂縫, 新斷層活動(dòng), 老斷層內(nèi)部斷裂縫分布基本不再變化, 斷層間具有較強(qiáng)的相互獨(dú)立性, 相互間的影響甚微。處于活動(dòng)中的斷層內(nèi), 新斷開(kāi)的粘結(jié)力鏈數(shù)、斷裂數(shù)多, 而停止活動(dòng)的斷層內(nèi)新斷開(kāi)的粘結(jié)力鏈則趨于零。該現(xiàn)象表明在變形過(guò)程中, 如果出現(xiàn)了新的斷層則老斷層基本停止活動(dòng)(孟令森等, 2007; 蔡申陽(yáng)等, 2016)。

    4.3 應(yīng)力和應(yīng)變

    收縮量不同時(shí)的累積體積應(yīng)變見(jiàn)圖10。由于左墻擠壓作用, 模型整體呈現(xiàn)體壓縮狀態(tài), 靠近擠壓端應(yīng)變較大, 遠(yuǎn)離擠壓端應(yīng)變較小, 斷層內(nèi)部體膨脹和體壓縮是共存的, 楔體淺部呈現(xiàn)體膨脹狀態(tài), 與前人研究結(jié)果一致(Morgan, 2015)。由于深部圍壓大, 模型體積膨脹受限, 深部體膨脹較少。

    圖7 每擠壓1 km過(guò)程中, 新斷開(kāi)的粘結(jié)力鏈數(shù)

    收縮量不同時(shí), 模型內(nèi)的累積變形應(yīng)變見(jiàn)圖11。變形應(yīng)變強(qiáng)的區(qū)域與體應(yīng)變較高的區(qū)域相互重疊(蔡申陽(yáng)等, 2016)。正向逆沖(紅色)往往伴隨有反向逆沖(藍(lán)色)斷層, 尤其當(dāng)?shù)撞炕搶訌?qiáng)度較弱時(shí), 這種共軛的“X”狀發(fā)育的斷層更為明顯(Morgan, 2015;李長(zhǎng)圣和尹宏偉, 2017; 李長(zhǎng)圣, 2019)。斷層F3(圖11c)和斷層F4(圖11d)的淺部變形應(yīng)變較大, 也表明斷層是從淺部開(kāi)始發(fā)育, 并逐步向深部發(fā)展。

    圖8 斷層的滑移曲線隨收縮量的變化(斷層F1、 F2、 F3和F4見(jiàn)圖4, 斷層滑移值選取綠色層測(cè)量)

    平均應(yīng)力演化過(guò)程見(jiàn)圖12, 平均應(yīng)力最大的地方均集中在靠近擠壓端的深部地層, 越遠(yuǎn)離擠壓端, 平均應(yīng)力越小, 地層越厚的地方平均應(yīng)力越大。收縮量為20 km時(shí), V1、V2、V3和V4處的平均應(yīng)力隨深度變化曲線見(jiàn)圖13。平均應(yīng)力與自重應(yīng)力分布基本一致, 地層厚度大的地方, 平均應(yīng)力大。選取水平方向H1測(cè)線, 可見(jiàn)平均應(yīng)力分布與地形起伏呈正相關(guān)關(guān)系(圖14)。

    隨著地層持續(xù)收縮, 最大剪切應(yīng)力首先在靠近擠壓端深部集中, 形成斷層F1, 之后不斷往擠壓前端傳播的, 持續(xù)在將要形成的新斷層F3處累積(圖15b), 直至斷層F3形成(圖15c), 剪切應(yīng)力開(kāi)始消散(圖15d), 繼續(xù)往前傳播, 在下一個(gè)將要形成的新斷層處累積。最大剪切應(yīng)力是一個(gè)體系的瞬時(shí)狀態(tài), 斷層形成后最大剪切應(yīng)力釋放, 斷層附近最大剪切應(yīng)力集中現(xiàn)象將消失。隨著擠壓進(jìn)行, 最大剪切應(yīng)力隨楔體逐漸向前傳播, 在推覆前端累積, 斷層形成, 應(yīng)力消散, 并在推覆前端繼續(xù)累積。盡管深部最大剪切應(yīng)力較大, 但是變形往往從淺部開(kāi)始(圖11),并逐步向深部擴(kuò)展。

    圖a中4個(gè)點(diǎn)分別選取自綠色標(biāo)志層中的斷層F1、F2、F3和F4內(nèi)部的4個(gè)顆粒, 為便于識(shí)別, 這4個(gè)顆粒半徑放大3倍, 藍(lán)色力鏈為監(jiān)測(cè)點(diǎn)周圍1 km范圍內(nèi)的粘結(jié)力鏈。

    藍(lán)色表示體壓縮, 紅色表示體膨脹, 顏色深淺表征體積變化的大小。其中, (e)中半透明黑色為標(biāo)志層。

    藍(lán)色表示逆時(shí)針剪切, 代表反沖斷層帶; 紅色表示順時(shí)針剪切, 代表逆沖斷層帶。顏色深淺表征了剪切應(yīng)變的大小。其中, (e)中黑色為標(biāo)志層。

    黑色線條為應(yīng)力等值線, 黑色點(diǎn)為|變形應(yīng)變|>4的點(diǎn), 表征了斷層的位置。子圖(e)中, V1、V2、V3和V4為應(yīng)力監(jiān)測(cè)線, 黃色點(diǎn)為PF1、PF2、PF3和PF4位置, 綠色層為標(biāo)志層。

    自重應(yīng)力按照ρgh計(jì)算, 其中, 巖層密度ρ取2500 km/m3, 重力加速度g取10, h為地層深度。監(jiān)測(cè)線V1、V2、V3、V4的位置見(jiàn)圖12, 黃色點(diǎn)為PF1、PF2、PF3和PF4的應(yīng)力值。

    5 結(jié) 論

    近年來(lái), 離散元方法越來(lái)越多地應(yīng)用到構(gòu)造模擬中。本文簡(jiǎn)述了離散元的基本原理, 通過(guò)雙軸實(shí)驗(yàn)標(biāo)定了顆粒的細(xì)觀參數(shù)對(duì)應(yīng)的地層宏觀力學(xué)參數(shù)(粘聚力10.5 MPa, 內(nèi)摩擦角18.6°), 實(shí)現(xiàn)了一個(gè)典型的擠壓構(gòu)造離散元數(shù)值模擬實(shí)驗(yàn), 并給出了基于離散元的構(gòu)造變形定量分析方法, 取得以下認(rèn)識(shí):

    (1) 該模擬實(shí)驗(yàn)的構(gòu)造以前展式的逆沖疊瓦斷層為主, 靠近擠壓端的斷層到遠(yuǎn)離擠壓端的斷層依次活動(dòng);

    (2) 裂縫與斷層形成有密切關(guān)系, 局部區(qū)域內(nèi)聚集的大量裂縫是產(chǎn)生斷層的誘因;

    (3) 斷層形成是瞬時(shí)的, 斷層形成初期, 粘結(jié)力鏈增量最大, 呈現(xiàn)開(kāi)口向上的弧形, 之后, 該新斷層將持續(xù)活動(dòng), 老斷層活動(dòng)基本停止, 新斷裂的粘結(jié)力鏈數(shù)逐漸減小; 斷層之間相互獨(dú)立, 斷層形成初期, 斷層內(nèi)物質(zhì)運(yùn)動(dòng)距離較小, 但是新產(chǎn)生裂縫數(shù)量最多; 斷層活動(dòng)后期, 斷層內(nèi)物質(zhì)隨著模型收縮, 位移呈現(xiàn)線形增加, 新產(chǎn)生的裂縫數(shù)較少。微斷裂往往從淺部發(fā)育并逐步向深部擴(kuò)展, 最終貫穿形成斷層;

    圖14 收縮量為20 km時(shí),水平方向平均應(yīng)力分布

    黑色線條為應(yīng)力等值線, 黑色點(diǎn)為|變形應(yīng)變|>4的點(diǎn), 表征了斷層的位置。圖(e)中, 綠色層為標(biāo)志層。

    (4) 體積應(yīng)變可以表征裂縫類型(拉張或壓縮), 變形應(yīng)變可以區(qū)分正向和反向逆沖斷層;

    (5) 隨著地層收縮, 平均應(yīng)力大小與地形起伏呈正相關(guān)關(guān)系, 呈現(xiàn)楔狀體分布, 靠近擠壓端地層厚度最大; 隨著地層持續(xù)收縮, 最大剪切應(yīng)力在遠(yuǎn)離擠壓端逐漸累積, 最大剪切應(yīng)力不斷往擠壓前端傳播的, 持續(xù)在將要形成的新斷層處累積, 直至該新斷層形成, 剪切應(yīng)力開(kāi)始消散, 繼續(xù)往前傳播, 在下一個(gè)將要形成的新斷層處累積。

    離散元方法將顆粒集合體模型視為若干離散單元的集合, 允許顆粒間產(chǎn)生大位移, 特別適合用于構(gòu)造變形定量研究。離散元方法在構(gòu)造變形、應(yīng)力應(yīng)變與裂縫預(yù)測(cè)定量研究中具有巨大潛力, 是未來(lái)的構(gòu)造變形定量研究的主要方向之一。

    致謝: 本文的數(shù)值計(jì)算在南京大學(xué)高性能計(jì)算中心的計(jì)算集群上完成, 數(shù)值模擬實(shí)驗(yàn)使用課題組自主研發(fā)的離散元數(shù)值模擬軟件完成, 更多構(gòu)造模擬示例見(jiàn)https: //geovbox.com。文中采用的應(yīng)變計(jì)算代碼, 修改自萊斯大學(xué)Julia K. Morgan和Thomas Fournier的腳本, 在此表示感謝。楔體要素測(cè)量程序源碼見(jiàn)https: //github.com/demsheng/Quantitative-wedge。中國(guó)石油大學(xué)(北京)余一欣副教授和劉志娜副教授對(duì)本文進(jìn)行了認(rèn)真而專業(yè)的審閱, 提出了建設(shè)性修改建議, 在此一并致以特別感謝。

    蔡申陽(yáng), 尹宏偉, 李長(zhǎng)圣, 賈東, 汪偉, 陳竹新, 魏東濤. 2016. 基于離散元數(shù)值模擬的應(yīng)變分析和裂縫預(yù)測(cè)技術(shù). 高校地質(zhì)學(xué)報(bào), 22(1): 183–193.

    李江海, 章雨, 王洪浩, 王殿舉. 2020. 庫(kù)車前陸沖斷帶西部古近系鹽構(gòu)造三維離散元數(shù)值模擬. 石油勘探與開(kāi)發(fā), 47(1): 65–76.

    李維波, 李江海, 王洪浩, 黃少英, 能源. 2017. 庫(kù)車前陸沖斷帶克拉蘇構(gòu)造帶變形影響因素分析——基于離散元數(shù)值模擬研究. 大地構(gòu)造與成礦學(xué), 41(6): 1001– 1010.

    李長(zhǎng)圣. 2014. 含礫滑帶土三軸剪切的精細(xì)數(shù)值模擬研究. 南京: 南京大學(xué)碩士學(xué)位論文.

    李長(zhǎng)圣. 2019. 基于離散元的褶皺沖斷帶構(gòu)造變形定量分析與模擬. 南京: 南京大學(xué)博士學(xué)位論文.

    李長(zhǎng)圣, 尹宏偉. 2017. 滑脫層強(qiáng)度對(duì)擠壓構(gòu)造的影響: 離散元數(shù)值模擬 // 2017中國(guó)地球科學(xué)聯(lián)合學(xué)術(shù)年會(huì)論文集. 北京: 中國(guó)和平音像電子出版社: 3879–3882.

    李長(zhǎng)圣, 尹宏偉, 劉春, 蔡申陽(yáng). 2017. 共享內(nèi)存式并行離散元程序的設(shè)計(jì)與測(cè)試. 南京大學(xué)學(xué)報(bào)(自然科學(xué)版), 53(6): 1161–1170.

    李長(zhǎng)圣, 張丹, 王宏憲, 獨(dú)莎莎. 2014. 基于CT掃描的土石混合體三維數(shù)值網(wǎng)格的建立. 巖土力學(xué), 35(9): 2731–2736.

    劉其鵬. 2010. 基于平均場(chǎng)理論的顆粒材料離散顆粒集合- Cosserat連續(xù)體模型多尺度模擬. 大連: 大連理工大學(xué)博士學(xué)位論文.

    劉源, Konietzky H. 2019. 基于離散元方法對(duì)走滑拉分盆地演化及次級(jí)斷裂擴(kuò)展過(guò)程研究. 地質(zhì)力學(xué)學(xué)報(bào), 25(5): 840–852.

    孟令森, 尹宏偉, 張潔, 徐士進(jìn). 2007. 巖石強(qiáng)度和應(yīng)變速率對(duì)水平擠壓變形影響的離散元模擬. 巖石學(xué)報(bào), 23(11): 2918–2926.

    沈禮, 賈東, 尹宏偉, 孫闖, 張勇, 范小根. 2012. 基于粒子成像測(cè)速(PIV)技術(shù)的褶皺沖斷構(gòu)造物理模擬. 地質(zhì)論評(píng), 58(3): 471–480.

    孫闖. 2017. 龍門山褶皺沖斷帶構(gòu)造物理模擬研究. 南京: 南京大學(xué)博士學(xué)位論文.

    吳珍云, 尹宏偉, 李長(zhǎng)圣, 孫業(yè)君, 李麗梅, 杜航, 何奕成, 劉博雅. 2019. 斷陷盆地正反轉(zhuǎn)構(gòu)造實(shí)驗(yàn)?zāi)M及對(duì)茅山斷裂帶的啟示. 南京大學(xué)學(xué)報(bào)(自然科學(xué)版), 55(5): 869–878.

    辛文, 陳漢林, 安凱旋, 張欲清, 楊樹(shù)鋒, 程曉敢, 林秀斌. 2020. 基于離散元數(shù)值模擬的西南天山山前沖斷帶構(gòu)造變形控制因素研究. 地質(zhì)學(xué)報(bào), 94(6): 1704–1715.

    徐雯嶠, 汪偉, 尹宏偉, 賈東, 李長(zhǎng)圣, 楊庚兄, 李剛. 2020. 庫(kù)車坳陷東西段鹽下構(gòu)造變形差異演化數(shù)值模擬分析. 地質(zhì)學(xué)報(bào), 94(6): 1740–1751.

    張潔, 尹宏偉, 孟令森, 徐士進(jìn). 2008. 主動(dòng)底辟鹽構(gòu)造的二維離散元模擬. 地球物理學(xué)進(jìn)展, 23(6): 1924– 1930.

    周易, 于福生, 劉志娜. 2019. 分層伸展疊加變形二維離散元模擬. 大地構(gòu)造與成礦學(xué), 43(2): 213–225.

    鄒瑋, 余一欣, 劉金水, 蔣一鳴, 唐賢君, 陳石, 余浪. 2021. 東海盆地西湖凹陷中央反轉(zhuǎn)構(gòu)造帶發(fā)育主控因素及寧波背斜形成過(guò)程. 石油學(xué)報(bào), 42(2): 176–185.

    Bigi S, Di Paolo L, Vadacca L, Gambardella G. 2010. Load and unload as interference factors on cyclical behavior and kinematics of Coulomb wedges: Insights from sandboxexperiments., 32(1): 28–44.

    Buiter S J. 2012. A review of brittle compressional wedge models., 530: 1–17.

    Buiter S J, Babeyko A Y, Ellis S, Gerya T V, Kaus B J, Kellner A, Schreurs G, Yamada Y. 2006. The numerical sandbox: Comparison of model results for a shortening and an extension experiment.,,, 253: 29–64.

    Buiter S J, Schreurs G, Albertz M, Gerya T V, Kaus B, Landry W, Le Pourhiet L, Mishin Y, Egholm D L, Cooke M. 2016. Benchmarking numerical models of brittle thrust wedges., 92: 140–177.

    Cui J, Jia D, Yin H W, Chen Z X, Li Y Q, Wang M M, Fan X G, Shen L, Sun C, Li Z G, Ma D L, Zhang Y K. 2020. The influence of a weak upper ductile detachment on the Longmen Shan fold-and-thrust belt (eastern margin of the Tibetan Plateau): Insights from sandbox experiments., 198, 104220.

    Cundall P A, Strack O D. 1979. A discrete numerical model for granular assemblies., 29(1): 47–65.

    Dahlen F A. 1984. Noncohesive critical Coulomb wedges: An exact solution.:, 89(B12): 10125–10133.

    Dahlen F A. 1990. Critical taper model of fold-and-thrust belts and accretionary wedges., 18(1): 55–99.

    Davis D, Suppe J, Dahlen F A. 1983. Mechanics of fold-and- thrust belts and accretionary wedges.:, 88(B2): 1153–1172.

    Dean S L, Morgan J K, Fournier T. 2013. Geometries of frontal fold and thrust belts: Insights from discrete element simulations., 53: 43–53.

    Gutscher M, Kukowski N, Malavieille J, Lallemand S. 1996. Cyclical behavior of thrust wedges: Insights from high basal friction sandbox experiments., 24(2): 135–138.

    Gutscher M, Kukowski N, Malavieille J, Lallemand S. 1998. Material transfer in accretionary wedges from analysis of a systematic series of analog experiments., 20(4): 407–416.

    Itasca Consulting Group. 2008. PFC2D (Particle Flow Code in 2 Dimensions) Online Manual (Version 4.0). Minnesota: Itasca Consulting Group Inc.

    Jaeger J C, Cook N G, Zimmerman R. 2009. Fundamentals of Rock Mechanics. John Wiley & Sons.

    Li C S, Yin H W, Jia D, Zhang J X, Wang W, Xu S J. 2018. Validation tests for discrete element codes using single- contact systems., 18(6), 6018011.

    Li C S, Yin H W, Wu C, Zhang Y C, Zhang J X, Wu Z Y, Wang W, Jia D, Guan S W, Ren R. 2021. Calibration of the discrete element method and modelling of shortening experiments., 9, 636512.

    Li C S, Zhang D, Du S S, Shi B. 2016. Computed tomography based numerical simulation for triaxial test of soil-rock mixture., 73: 179–188.

    Liu Z, Koyi H A, Swantesson J O, Nilfouroushan F, Reshetyuk Y. 2013. Kinematics and 3-D internal deformation of granular slopes: Analogue models and natural landslides., 53: 27–42.

    Malvern L E. 1969. Introduction to the Mechanics of a Continuous Medium. New Jersey, United States. Prentice- Hall, Inc.

    Mathworks. 2015. MATLAB documentation. Natick, United States. The MathWorks Inc.

    Means W D. 2012. Stress and strain: Basic concepts of continuum mechanics for geologists. Springer Science & Business Media.

    Morgan J K. 2015. Effects of cohesion on the structural and mechanical evolution of fold and thrust belts and contractional wedges: Discrete element simulations.:, 120(5): 3870– 3896.

    Oertel G. 1996. Stress and deformation: A handbook on tensors in geology. Oxford University Press.

    O’Sullivan C. 2011. Particulate Discrete Element Modelling: A Geomechanics Perspective. Taylor & Francis.

    Schreurs G, Buiter S J, Boutelier D, Corti G, Costa E, Cruden A R, Daniel J, Hoth S, Koyi H A, Kukowski N. 2006. Analogue benchmarks of shortening and extension experiments.,,, 253: 1–27.

    Schreurs G, Buiter S J, Boutelier J, Burberry C, Callot J, Cavozzi C, Cerca M, Chen J, Cristallini E, Cruden A R. 2016. Benchmarking analogue models of brittle thrust wedges., 92: 116–139.

    Sun C, Jia D, Yin H W, Chen Z X, Li Z G, Shen L, Wei D T, Li Y Q, Yan B, Wang M M, Fang S Z, Cui J. 2016. Sandbox modeling of evolving thrust wedges with different preexisting topographic relief: Implications for the Longmen Shan thrust belt, eastern Tibet.:, 121(6): 4591–4614.

    Suppe J. 2007. Absolute fault and crustal strength from wedge tapers., 35(12): 1127–1130.

    Wang X Y, Morgan J, Bangs N. 2021. Relationships among forearc structure, fault slip, and earthquake magnitude: Numerical simulations with applications to the Central Chilean Margin., e2021GL092521. doi:10.1002/essoar.10506057.2

    Xu W Q, Yin H W, Jia D, Li C S, Wang W, Yang G X, He W H, Chen Z X, Ren R. 2021. Structural features and evolution of the northwestern Sichuan Basin: Insights from discrete numerical simulations., 9, 653395.

    Yin H W, Zhang J, Meng L S, Liu Y P, Xu S. 2009. Discrete element modeling of the faulting in the sedimentary cover above an active salt diapir., 31(9): 989–995.

    Quantitative Analysis and Simulation of Compressive Tectonics Based on Discrete Element Method

    LI Changsheng1, 2, 3, 4, YIN Hongwei3*, XU Wenqiao3, WU Zhenyun1, 2, GUAN Shuwei4, JIA Dong3, REN Rong4

    (1. State Key Laboratory of Nuclear Resources and Environment, Nanchang 330013, Jiangxi, China; 2. School of Earth Sciences, East China University of Technology, Nanchang 330013, Jiangxi, China; 3. School of Earth Sciences and Engineering, Nanjing University, Nanjing 210023, Jiangsu, China; 4. PetroChina Research Institute of Petroleum Exploration & Development, Beijing 100083, China)

    With development of discrete element theory and computer technology, the discrete element method has been widely used in structural simulation at different scales. Compared with the traditional scaled analog (sandbox) experiments, the discrete element numerical simulation can precisely control the experimental boundary conditions and quantitatively analyze structural deformation process, which is helpful for deeply understanding the mechanical properties and deformation mechanism of stratum from the mesoscopic scale. In this paper, based on a typical discrete element simulation of compressive tectonics, we systematically expounded the tectonic deformation quantitative analysis method, highlighted the superiority of tectonic deformation in mesoscale, and obtained the following results: (1) In the typical numerical simulation of compressive tectonics, the forward spreading thrust imbricated faults are dominant, and the faults generate in sequence from the compression end to those far from the compression end. (2) Fractures are closely related to the formation of faults, the interval in which the fracture generation reaches its peak precedes the one in which the fault appears the most active. (3) In the early stage of fault formation, the faults have a small displacement. However, the number of new fractures is the largest in this stage. In the later stage of fault activity, the fault-slip extension amount increases linearly with the shrinkage of the model, and the number of new fractures is less. (4) Volumetric strain can represent fracture types (tension or compression), and deformation strain can distinguish forward and back thrust faults. (5) The average stress is positively correlated with topographic relief, and the maximum shear stress continues to accumulate under the new fault to be formed until the new fault is formed, then the shear stress begins to dissipate and spread forward, accumulating at the next new fault to be formed. The results of this study show that the discrete element method has great potential in the quantitative research of structural deformation, stress-strain, and fracture prediction.

    discrete element method; structural deformation; compressive tectonics; stress and strain; quantitative analysis

    P542; TE352

    A

    1001-1552(2022)04-0645-017

    10.16539/j.ddgzyckx.2022.04.001

    2020-12-06;

    2021-06-16

    國(guó)家自然科學(xué)基金項(xiàng)目(42102270、41972219、42172232、41927802、41572187、41602208)、東華理工大學(xué)博士啟動(dòng)基金項(xiàng)目(DHBK2019024、DHBK2019053)和中國(guó)石油天然氣股份有限公司項(xiàng)目(2018A-0101)聯(lián)合資助。

    李長(zhǎng)圣(1989–), 男, 博士, 主要從事構(gòu)造及巖土體相關(guān)數(shù)值模擬的研究工作。E-mail: sheng0619@163.com

    尹宏偉(1971–), 男, 教授, 主要從事構(gòu)造分析和構(gòu)造物理、數(shù)值模擬等領(lǐng)域的研究工作。E-mail: hwyin@nju.edu.cn

    猜你喜歡
    變形模型
    一半模型
    重要模型『一線三等角』
    談詩(shī)的變形
    重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
    “我”的變形計(jì)
    變形巧算
    例談拼圖與整式變形
    會(huì)變形的餅
    3D打印中的模型分割與打包
    FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
    麻豆国产av国片精品| 亚洲第一区二区三区不卡| 欧美变态另类bdsm刘玥| 九九爱精品视频在线观看| 亚洲成a人片在线一区二区| 午夜精品在线福利| 一个人观看的视频www高清免费观看| 噜噜噜噜噜久久久久久91| 爱豆传媒免费全集在线观看| 最好的美女福利视频网| 日本-黄色视频高清免费观看| 亚洲经典国产精华液单| 久久久久久久久久成人| av专区在线播放| 中文字幕人妻熟人妻熟丝袜美| 成人鲁丝片一二三区免费| 99在线视频只有这里精品首页| 1024手机看黄色片| 亚洲精品影视一区二区三区av| 一区二区三区四区激情视频 | 亚洲av电影不卡..在线观看| 日韩一区二区三区影片| 亚洲国产色片| 免费观看的影片在线观看| 又粗又硬又长又爽又黄的视频 | 国产老妇伦熟女老妇高清| 国内精品美女久久久久久| 亚洲精品久久久久久婷婷小说 | 久久久久性生活片| 亚洲人与动物交配视频| 男人的好看免费观看在线视频| 亚洲精品影视一区二区三区av| 久久午夜福利片| 岛国在线免费视频观看| 免费大片18禁| 一级二级三级毛片免费看| 在线观看美女被高潮喷水网站| 成年女人永久免费观看视频| 国产男人的电影天堂91| avwww免费| 国产一区二区在线观看日韩| 久久国产乱子免费精品| 蜜臀久久99精品久久宅男| 99国产极品粉嫩在线观看| 欧美又色又爽又黄视频| 国产亚洲精品久久久久久毛片| 亚洲中文字幕一区二区三区有码在线看| 国产精品乱码一区二三区的特点| 国产精品麻豆人妻色哟哟久久 | 欧美bdsm另类| 免费观看a级毛片全部| 中文欧美无线码| 午夜福利在线在线| 99国产极品粉嫩在线观看| 日产精品乱码卡一卡2卡三| 日本免费一区二区三区高清不卡| 婷婷精品国产亚洲av| 熟女电影av网| 国产免费一级a男人的天堂| 成人毛片60女人毛片免费| 日韩av不卡免费在线播放| 人妻夜夜爽99麻豆av| 两个人的视频大全免费| 寂寞人妻少妇视频99o| 亚洲欧美精品自产自拍| 欧美最新免费一区二区三区| 变态另类成人亚洲欧美熟女| 大香蕉久久网| 网址你懂的国产日韩在线| 看十八女毛片水多多多| 国产毛片a区久久久久| 亚洲中文字幕日韩| 不卡一级毛片| 国产乱人视频| 久久精品国产自在天天线| 久久精品国产清高在天天线| 国产伦精品一区二区三区视频9| www.色视频.com| 精品午夜福利在线看| 国产高清视频在线观看网站| 麻豆国产97在线/欧美| 国产一区二区亚洲精品在线观看| 少妇裸体淫交视频免费看高清| 久久久久久大精品| 我要看日韩黄色一级片| 中文精品一卡2卡3卡4更新| 中文字幕久久专区| 中文字幕免费在线视频6| 午夜精品一区二区三区免费看| 熟妇人妻久久中文字幕3abv| 中文亚洲av片在线观看爽| av在线观看视频网站免费| 久久久久久伊人网av| 亚洲av一区综合| 国产亚洲精品久久久com| 久久精品国产99精品国产亚洲性色| 亚洲欧美日韩高清在线视频| www.av在线官网国产| 边亲边吃奶的免费视频| 久久精品久久久久久噜噜老黄 | 免费观看人在逋| 国产精华一区二区三区| 日本三级黄在线观看| 亚洲欧美成人精品一区二区| 亚洲乱码一区二区免费版| 日韩成人av中文字幕在线观看| 91精品国产九色| h日本视频在线播放| 亚洲综合色惰| 午夜精品在线福利| 少妇的逼水好多| 久久久午夜欧美精品| 欧美性猛交╳xxx乱大交人| 美女高潮的动态| 少妇熟女aⅴ在线视频| 一区福利在线观看| 在线观看午夜福利视频| 久久久久久久亚洲中文字幕| 亚洲激情五月婷婷啪啪| 嫩草影院新地址| 国产91av在线免费观看| .国产精品久久| 国产精品人妻久久久影院| 国产精品久久久久久精品电影小说 | 一级黄片播放器| 国产成人午夜福利电影在线观看| 村上凉子中文字幕在线| 欧美激情在线99| 午夜久久久久精精品| 亚洲欧美精品自产自拍| 亚洲最大成人av| 特大巨黑吊av在线直播| 变态另类成人亚洲欧美熟女| 一边亲一边摸免费视频| 少妇熟女欧美另类| 两性午夜刺激爽爽歪歪视频在线观看| 久久亚洲精品不卡| 久久国内精品自在自线图片| 人人妻人人看人人澡| 亚洲成人中文字幕在线播放| 99riav亚洲国产免费| 免费观看精品视频网站| 国产av在哪里看| 亚洲aⅴ乱码一区二区在线播放| 久久久久久久久久久丰满| 午夜福利在线观看免费完整高清在 | 中文字幕免费在线视频6| 熟女电影av网| 亚洲无线在线观看| 成人亚洲精品av一区二区| 亚洲中文字幕日韩| 淫秽高清视频在线观看| 日本黄色片子视频| 亚洲国产精品合色在线| 我的女老师完整版在线观看| 午夜爱爱视频在线播放| 欧美人与善性xxx| 国产伦精品一区二区三区视频9| 亚洲欧美日韩高清专用| 欧美日韩国产亚洲二区| 亚洲欧美成人综合另类久久久 | av在线老鸭窝| 国产亚洲5aaaaa淫片| 欧美日韩综合久久久久久| 身体一侧抽搐| 久久久精品欧美日韩精品| 精品日产1卡2卡| 国语自产精品视频在线第100页| 亚洲国产欧洲综合997久久,| 青春草国产在线视频 | 欧美性感艳星| 亚洲天堂国产精品一区在线| 男女啪啪激烈高潮av片| 亚洲av熟女| 最近最新中文字幕大全电影3| 1000部很黄的大片| 免费大片18禁| 老师上课跳d突然被开到最大视频| 国产乱人视频| 嫩草影院入口| 免费观看精品视频网站| 18禁黄网站禁片免费观看直播| 日韩欧美精品免费久久| 看片在线看免费视频| 免费观看a级毛片全部| 久久九九热精品免费| 国产伦理片在线播放av一区 | 99久国产av精品| 九九在线视频观看精品| 色哟哟·www| 中文字幕制服av| 日本-黄色视频高清免费观看| 18禁在线播放成人免费| 精品久久久久久久久久久久久| 日韩人妻高清精品专区| 国产精品久久久久久av不卡| 国产精品电影一区二区三区| 26uuu在线亚洲综合色| 久久九九热精品免费| 亚洲,欧美,日韩| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 国产精品一区二区在线观看99 | 中文字幕av成人在线电影| 日本av手机在线免费观看| 99久久成人亚洲精品观看| 波多野结衣高清无吗| av在线天堂中文字幕| 九九爱精品视频在线观看| 久久热精品热| 美女xxoo啪啪120秒动态图| 国产精品一二三区在线看| 成人高潮视频无遮挡免费网站| 国内少妇人妻偷人精品xxx网站| 一级毛片电影观看 | 国产色婷婷99| 亚洲精品日韩av片在线观看| 国内精品宾馆在线| 欧美人与善性xxx| av又黄又爽大尺度在线免费看 | 亚洲人成网站在线播| 国产爱豆传媒在线观看| 六月丁香七月| 国产精品不卡视频一区二区| 久久久久久久久大av| 日本撒尿小便嘘嘘汇集6| 熟女电影av网| 一个人免费在线观看电影| 亚洲欧美日韩高清在线视频| 午夜福利在线在线| 99久久精品国产国产毛片| 国产免费男女视频| 欧美精品国产亚洲| 国产一级毛片在线| 国产女主播在线喷水免费视频网站 | 欧美日本亚洲视频在线播放| 美女cb高潮喷水在线观看| 国产伦在线观看视频一区| 国产精品人妻久久久影院| 女的被弄到高潮叫床怎么办| 日本免费一区二区三区高清不卡| videossex国产| 91狼人影院| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 成年女人看的毛片在线观看| 亚洲av第一区精品v没综合| 成人一区二区视频在线观看| 成人国产麻豆网| 精品无人区乱码1区二区| 亚洲中文字幕日韩| 欧美激情在线99| 午夜免费男女啪啪视频观看| 久久精品国产清高在天天线| 极品教师在线视频| 在线观看av片永久免费下载| 中国美女看黄片| 中文资源天堂在线| 五月玫瑰六月丁香| 国产精品永久免费网站| 美女大奶头视频| 精品欧美国产一区二区三| 看十八女毛片水多多多| 久久综合国产亚洲精品| 国产色婷婷99| 免费av毛片视频| 免费看日本二区| 国内少妇人妻偷人精品xxx网站| videossex国产| 欧美激情国产日韩精品一区| 久久久久久久午夜电影| 国产大屁股一区二区在线视频| 精华霜和精华液先用哪个| 亚洲人与动物交配视频| 成年女人看的毛片在线观看| 亚洲精品影视一区二区三区av| 成人永久免费在线观看视频| 国产成人福利小说| 国产精品麻豆人妻色哟哟久久 | 欧美日韩在线观看h| 精品人妻熟女av久视频| 午夜福利在线观看免费完整高清在 | 国产午夜精品一二区理论片| 国产精品一区二区三区四区久久| 九九爱精品视频在线观看| 久久人人精品亚洲av| 久99久视频精品免费| 精品不卡国产一区二区三区| 亚洲aⅴ乱码一区二区在线播放| 国产亚洲精品av在线| 五月玫瑰六月丁香| 日韩欧美国产在线观看| 久久精品夜夜夜夜夜久久蜜豆| 日本黄色片子视频| 国产三级中文精品| 亚洲国产精品成人久久小说 | 国产成人aa在线观看| 国产老妇女一区| 亚洲aⅴ乱码一区二区在线播放| 久久亚洲国产成人精品v| 亚洲内射少妇av| 真实男女啪啪啪动态图| 久久精品国产亚洲av香蕉五月| 美女高潮的动态| 亚洲激情五月婷婷啪啪| 亚洲欧美中文字幕日韩二区| 夜夜看夜夜爽夜夜摸| 天天躁日日操中文字幕| 不卡视频在线观看欧美| 中文字幕熟女人妻在线| 一级毛片aaaaaa免费看小| 99热网站在线观看| 久久精品国产亚洲av天美| 日本与韩国留学比较| 亚洲精品久久国产高清桃花| 国产精品国产高清国产av| 黄色视频,在线免费观看| 免费一级毛片在线播放高清视频| 亚州av有码| 又爽又黄a免费视频| 婷婷色综合大香蕉| 嘟嘟电影网在线观看| 亚洲激情五月婷婷啪啪| 国产精品一区www在线观看| 中文亚洲av片在线观看爽| 欧美成人a在线观看| 亚洲国产精品合色在线| 校园人妻丝袜中文字幕| 欧美激情国产日韩精品一区| 欧美一区二区亚洲| 色5月婷婷丁香| 国产精品日韩av在线免费观看| 男女那种视频在线观看| 在线免费观看的www视频| 99久久成人亚洲精品观看| 久久精品夜夜夜夜夜久久蜜豆| 成人二区视频| 亚洲婷婷狠狠爱综合网| 看免费成人av毛片| 午夜亚洲福利在线播放| 久久亚洲国产成人精品v| 国产精品爽爽va在线观看网站| 国产精品电影一区二区三区| 一边亲一边摸免费视频| 亚洲国产精品久久男人天堂| 美女大奶头视频| 久久久国产成人精品二区| 波多野结衣高清无吗| 欧美xxxx性猛交bbbb| 成人毛片a级毛片在线播放| 国产精华一区二区三区| 久久精品91蜜桃| 好男人视频免费观看在线| 黄色一级大片看看| 两个人视频免费观看高清| 国产精品99久久久久久久久| 18禁在线播放成人免费| 国产成人一区二区在线| 免费观看人在逋| 伦精品一区二区三区| 日韩强制内射视频| 人人妻人人看人人澡| 麻豆国产97在线/欧美| 免费看日本二区| 精品人妻视频免费看| 99热这里只有是精品在线观看| 蜜臀久久99精品久久宅男| 久久精品国产亚洲网站| 久久99蜜桃精品久久| 国产69精品久久久久777片| 波多野结衣高清无吗| 99热6这里只有精品| 久久久午夜欧美精品| 少妇熟女欧美另类| 国内少妇人妻偷人精品xxx网站| 亚洲最大成人手机在线| 亚洲中文字幕一区二区三区有码在线看| 国产av在哪里看| 国产视频内射| 亚洲自拍偷在线| 乱系列少妇在线播放| 精品人妻偷拍中文字幕| 美女 人体艺术 gogo| 精品人妻偷拍中文字幕| 久久精品国产亚洲网站| 91aial.com中文字幕在线观看| 青春草国产在线视频 | 1024手机看黄色片| 国内少妇人妻偷人精品xxx网站| 国产美女午夜福利| av在线观看视频网站免费| 美女高潮的动态| 亚洲国产精品sss在线观看| 国产v大片淫在线免费观看| eeuss影院久久| 成熟少妇高潮喷水视频| 日韩中字成人| 亚洲人成网站在线播| 国产蜜桃级精品一区二区三区| 天堂影院成人在线观看| 波多野结衣高清作品| 亚洲丝袜综合中文字幕| 嫩草影院精品99| 国产精品不卡视频一区二区| 日本爱情动作片www.在线观看| 精品少妇黑人巨大在线播放 | 麻豆av噜噜一区二区三区| 国产真实伦视频高清在线观看| 女人被狂操c到高潮| 国产 一区精品| 极品教师在线视频| 精品一区二区免费观看| 久久久精品大字幕| 国产一区二区在线av高清观看| 免费观看精品视频网站| a级毛片a级免费在线| 欧美高清成人免费视频www| 噜噜噜噜噜久久久久久91| 欧美人与善性xxx| 欧美激情在线99| 免费av不卡在线播放| 男人的好看免费观看在线视频| 身体一侧抽搐| 少妇熟女欧美另类| eeuss影院久久| 国内久久婷婷六月综合欲色啪| 色吧在线观看| 国产伦精品一区二区三区四那| or卡值多少钱| 九草在线视频观看| 一本久久中文字幕| 夫妻性生交免费视频一级片| 一本久久中文字幕| 日本免费一区二区三区高清不卡| 啦啦啦啦在线视频资源| 国产不卡一卡二| 别揉我奶头 嗯啊视频| 爱豆传媒免费全集在线观看| 国产色爽女视频免费观看| 3wmmmm亚洲av在线观看| 尤物成人国产欧美一区二区三区| 国产男人的电影天堂91| 中国美白少妇内射xxxbb| 国产黄片视频在线免费观看| .国产精品久久| 99久国产av精品国产电影| 男的添女的下面高潮视频| a级毛色黄片| 欧美日本亚洲视频在线播放| 久久精品国产亚洲av天美| 日韩三级伦理在线观看| 日韩欧美在线乱码| 亚洲国产精品成人综合色| 欧美日韩综合久久久久久| 国产片特级美女逼逼视频| 亚洲精品乱码久久久久久按摩| 桃色一区二区三区在线观看| 久久这里只有精品中国| 99热精品在线国产| 欧美性猛交╳xxx乱大交人| 国产精品国产高清国产av| 高清毛片免费观看视频网站| 日韩av不卡免费在线播放| 久久精品国产99精品国产亚洲性色| 久久这里只有精品中国| 联通29元200g的流量卡| 午夜亚洲福利在线播放| 国产亚洲精品久久久久久毛片| 国产人妻一区二区三区在| 蜜臀久久99精品久久宅男| 午夜精品国产一区二区电影 | 欧美高清性xxxxhd video| 99热只有精品国产| 午夜视频国产福利| 亚洲av第一区精品v没综合| 在线观看免费视频日本深夜| 亚洲欧美中文字幕日韩二区| 免费观看a级毛片全部| 一本精品99久久精品77| 九九热线精品视视频播放| 黄色配什么色好看| 日本一二三区视频观看| 国产av在哪里看| 内地一区二区视频在线| 床上黄色一级片| 精品少妇黑人巨大在线播放 | 国产色爽女视频免费观看| 久久久久久久午夜电影| 99久久九九国产精品国产免费| 久久九九热精品免费| 亚洲真实伦在线观看| 欧美另类亚洲清纯唯美| 国产成人精品婷婷| 一个人看视频在线观看www免费| 日韩高清综合在线| 大型黄色视频在线免费观看| 嘟嘟电影网在线观看| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 久久99热这里只有精品18| 免费人成在线观看视频色| 日韩成人av中文字幕在线观看| 美女脱内裤让男人舔精品视频 | 国内精品美女久久久久久| 亚洲成a人片在线一区二区| 99久久精品热视频| 国产午夜精品一二区理论片| 性插视频无遮挡在线免费观看| 日本黄大片高清| 国产视频内射| av免费在线看不卡| 日日摸夜夜添夜夜添av毛片| 亚洲婷婷狠狠爱综合网| 久久精品人妻少妇| 成年av动漫网址| 床上黄色一级片| 看免费成人av毛片| 久久久国产成人免费| 亚洲国产精品sss在线观看| 久久久a久久爽久久v久久| 听说在线观看完整版免费高清| 国产白丝娇喘喷水9色精品| 老司机影院成人| 欧美日本亚洲视频在线播放| 最近最新中文字幕大全电影3| 嫩草影院精品99| 国产高清不卡午夜福利| 久久精品国产亚洲网站| 亚洲av男天堂| 国产伦理片在线播放av一区 | 国产精品国产高清国产av| 1024手机看黄色片| 3wmmmm亚洲av在线观看| 亚洲丝袜综合中文字幕| 少妇猛男粗大的猛烈进出视频 | 亚洲欧美精品综合久久99| 免费电影在线观看免费观看| av在线亚洲专区| 春色校园在线视频观看| 日韩欧美精品v在线| 亚洲av成人精品一区久久| 一级毛片久久久久久久久女| 国产亚洲精品久久久久久毛片| 99热这里只有是精品在线观看| 高清日韩中文字幕在线| 欧美日韩国产亚洲二区| 久久九九热精品免费| 亚洲国产精品sss在线观看| 久久久精品大字幕| 两个人视频免费观看高清| 国产高潮美女av| 欧美成人精品欧美一级黄| 99在线人妻在线中文字幕| 婷婷亚洲欧美| 三级毛片av免费| 国产高清激情床上av| 国内精品久久久久精免费| 听说在线观看完整版免费高清| 熟女电影av网| 亚洲精品456在线播放app| 欧美精品国产亚洲| 欧美不卡视频在线免费观看| 久久精品夜色国产| 国产精品一及| 在线播放无遮挡| 禁无遮挡网站| 亚洲乱码一区二区免费版| 国产日韩欧美在线精品| 欧美zozozo另类| 成人毛片a级毛片在线播放| 麻豆成人午夜福利视频| 日日啪夜夜撸| 高清在线视频一区二区三区 | 亚洲精品成人久久久久久| 久久久久久大精品| 99视频精品全部免费 在线| 国产日本99.免费观看| 亚洲av男天堂| 精品少妇黑人巨大在线播放 | 国产激情偷乱视频一区二区| 一区二区三区免费毛片| 黄片无遮挡物在线观看| 大又大粗又爽又黄少妇毛片口| 久久韩国三级中文字幕| 亚洲成人av在线免费| 狂野欧美白嫩少妇大欣赏| 我的女老师完整版在线观看| 天天一区二区日本电影三级| 最近2019中文字幕mv第一页| 91麻豆精品激情在线观看国产| 国产免费男女视频| 精品人妻熟女av久视频| 亚洲av男天堂| 伊人久久精品亚洲午夜| 国产精品久久久久久精品电影小说 | 色吧在线观看| 日本熟妇午夜| 51国产日韩欧美| 美女国产视频在线观看| av在线亚洲专区| 亚洲国产精品成人久久小说 | 欧美变态另类bdsm刘玥| 18禁裸乳无遮挡免费网站照片| 97人妻精品一区二区三区麻豆| 国产激情偷乱视频一区二区| 亚洲精品国产av成人精品| 国产一区二区在线观看日韩| 在线国产一区二区在线| 一本久久精品| 久久精品国产自在天天线| 欧美高清成人免费视频www| 99久国产av精品国产电影| 亚洲欧洲国产日韩| 亚洲人成网站高清观看| 日韩成人av中文字幕在线观看| 大型黄色视频在线免费观看| 国产色爽女视频免费观看| 亚洲国产精品成人综合色|