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

    凍融循環(huán)作用下露天煤礦內(nèi)排土場土體力學特征及強度劣化機理

    2023-12-05 05:43:48張合勇王雪冬朱永東王海鵬漆利輝
    煤田地質(zhì)與勘探 2023年11期
    關(guān)鍵詞:模型

    張合勇,王雪冬,朱永東,王海鵬,漆利輝

    (遼寧工程技術(shù)大學 礦業(yè)學院,遼寧 阜新 123000)

    我國露天煤礦的采礦規(guī)模持續(xù)擴大,開采過程中產(chǎn)生的棄土礦渣形成了大量平臺緊實、邊坡松散高陡的排土場[1]。由于我國大部分露天煤礦位于高寒生態(tài)脆弱地區(qū),受凍融作用影響,土體的力學性質(zhì)會隨著溫度和含水率的變化表現(xiàn)出復雜的動態(tài)特性[2],常導致排土場凍脹、融凍滑塌、剝蝕等變形破壞。因此,研究凍融循環(huán)作用下土體力學特征及強度劣化機理對礦山安全生產(chǎn)和生態(tài)修復等工作具有十分重要的價值。

    近年來,凍融循環(huán)作用影響土體物理力學性質(zhì)的研究主要集中在不同凍融循環(huán)次數(shù)中應(yīng)力-變形、抗剪強度、黏聚力和內(nèi)摩擦角等方面的變化規(guī)律[3]。隨著凍融次數(shù)的增加,土體的抗剪強度和黏聚力多數(shù)表現(xiàn)出先降低后保持穩(wěn)定或先降低后增加而后保持穩(wěn)定的規(guī)律[4-5],內(nèi)摩擦角的變化則沒有明顯的規(guī)律性;應(yīng)力-應(yīng)變曲線表現(xiàn)出不受影響或由軟化型向硬化型過渡的變化[6-7]。為了研究凍融循環(huán)對土體內(nèi)部的影響,采用核磁共振測試、掃描電鏡觀測和壓汞試驗等[8-9],從微觀角度觀察到土體的變化不只有土顆粒的凍縮和水的相變[10],還存在水分運移。在水分場、溫度場的相互作用及變化規(guī)律的控制下,土體原本的穩(wěn)定結(jié)構(gòu)被破壞,造成強度劣化。

    國內(nèi)外學者通過對凍土內(nèi)部水分場、溫度場和應(yīng)力場的研究,分析引起土體強度劣化的機理,建立了多種條件下的理論模型[11-13],但理論模型難以說明水、熱、力三場相互作用下土體的劣化過程,而數(shù)值模擬則可以彌補理論模型的不足[14]??紤]到離散元法可以通過顆粒流理論和有限差分法等思想較好地從細觀角度模擬土的力學行為和內(nèi)部結(jié)構(gòu)的相關(guān)變化,桑宏偉等[15]采用離散元法實現(xiàn)熱量在顆粒間的傳遞,模擬出能源樁與周圍土層換熱的過程。尹楠等[16]基于離散元法模擬不同圍壓、不同溫度下凍土的靜三軸試驗。Le Tiancheng 等[17]基于離散元軟件,模擬出土樣蒸發(fā)、收縮和開裂的過程,均取得了良好的研究成果。因此,筆者選擇離散元法模擬水分場、溫度場和應(yīng)力場的相互作用來分析凍融循環(huán)作用下土體強度的劣化機理。

    筆者以內(nèi)蒙古元寶山露天煤礦內(nèi)排土場典型黏土為研究對象,通過進行天然含水率下不同凍融循環(huán)次數(shù)的室內(nèi)試驗,獲取凍融作用對土料宏觀力學性質(zhì)的影響規(guī)律。采用離散元數(shù)值模擬,研究凍融循環(huán)中溫度場、水分場的變化導致土體強度劣化的影響機理,以期為凍區(qū)排土場邊坡穩(wěn)定性分析和其他工程建設(shè)提供一定的參考。

    1 試驗材料與試驗方案

    1.1 試驗材料

    試樣所用土料取自內(nèi)蒙古元寶山露天煤礦東南幫內(nèi)排土場,取土地理位置及顆粒級配曲線如圖1 所示。按照試驗規(guī)程[18]進行相關(guān)室內(nèi)土性試驗得到土料的天然干密度為1.63 g/cm3,天然含水率為17.6%,液限和塑限分別為33.8%、19.3%,塑性指數(shù)為14.5。通過塑性指數(shù)判定土料類型為黏土,根據(jù)顆粒級配、細顆粒含量和液塑限認為該土料屬于凍融敏感性土[19]。

    圖1 取土位置和顆粒級配曲線Fig.1 Sampling location and particle gradation curve

    1.2 試樣制備

    由于排土場是人工堆積而成的松散體,具有質(zhì)地疏松、強度較低等特點。土料的離散性較大,為獲取具有代表性的土樣,在多個平臺上收集一定量的原狀土料,通過室內(nèi)試驗測得土料的天然含水率和天然干密度,依據(jù)上述條件進行重塑,制備直徑39.1 mm、高度80 mm 的三軸試樣。為了探究無外界水分交換條件下凍融作用引起的土體內(nèi)部水分運移及相變對力學性質(zhì)的影響,試驗采用保鮮膜將試樣包裹后再放入密封袋,模擬無外界水分交換的封閉系統(tǒng)。具體重塑和制樣步驟及過程如圖2 所示。

    圖2 土料重塑、制樣步驟及圖片F(xiàn)ig.2 Soil remolded and sample preparation steps and pictures

    1.3 試驗方案

    試驗采用的凍融設(shè)備包括冰箱和恒溫恒濕箱。冰箱最低凍結(jié)溫度為-40℃;恒溫恒濕箱溫控范圍0~50℃,精度±0.1℃??紤]到研究區(qū)2011-2022 年中最低的月平均氣溫為-18℃,本試驗條件設(shè)置為:凍結(jié)溫度-20℃、融化溫度+20℃,凍結(jié)和融化時間均為12 h,即24 h 為一個周期。試驗共設(shè)置7 組,每組設(shè)置3 個試樣及1 個備用,共需28 個,循環(huán)次數(shù)分別為0、1、3、5、10、15、25 次。試樣達到凍融次數(shù)后采用TSZ10-1 型三軸儀進行不固結(jié)不排水試驗,試驗設(shè)置圍壓分別為100、200 和300 kPa,軸向應(yīng)變速率設(shè)為0.8 mm/min,當應(yīng)變量達到15%時停止剪切,取峰值為最大破壞強度。

    2 凍融循環(huán)試驗結(jié)果與分析

    2.1 對試樣質(zhì)量和宏觀變形的影響

    試樣在凍融過程中質(zhì)量的變化從圖3 可以看出,試樣質(zhì)量的損失率很低,雖然呈緩慢增長,但最終損失率僅為0.124%,即損失0.23 g 水分,說明在試驗中采用的密封處理方法有效避免了水分流失對試驗結(jié)果的影響。通過高度和直徑變化可以看出,首次凍融后兩者的凍縮均相對明顯,分別為0.095%和0.347%;高度的凍縮率在前期增長較快,至第5 次時增加到0.276%,隨后變化減緩并基本穩(wěn)定在0.3%;而直徑凍縮率經(jīng)過首次顯著變化后,開始緩慢增長并逐漸趨于穩(wěn)定。經(jīng)分析,在凍融循環(huán)試驗中,凍結(jié)時土骨架收縮,因土體含水率較低,此過程以孔隙體積被壓縮為主;融化后土骨架的收縮無法完全恢復,所以宏觀上表現(xiàn)為凍縮,此外,土在凍結(jié)過程中形成的冷生構(gòu)造會產(chǎn)生微裂隙[20],且試樣在試驗時為平放,受到重力影響,所以直徑收縮率大于高度收縮率;而高度方向的變化受重力影響較小,主要是凍融作用的結(jié)果,說明該土體受凍融作用的影響主要體現(xiàn)在前期。首次體積凍縮率最大,達到0.786%,隨著凍融次數(shù)增加,增長逐漸減緩并趨于穩(wěn)定;可以看出,由于土體的凍融敏感性,首次凍融作用下土體結(jié)構(gòu)、顆粒間的膠結(jié)和孔隙變化幅度最大[21],造成的體積縮小最為明顯,隨著凍融次數(shù)增加,凍融作用的影響減小,內(nèi)部結(jié)構(gòu)逐漸調(diào)整完成,體積趨于穩(wěn)定。

    圖3 試樣的質(zhì)量和宏觀變形與凍融次數(shù)的關(guān)系Fig.3 The relationship between the quality and macroscopic deformation of the sample and the number of freeze-thaw cycles

    2.2 對土體力學性質(zhì)的影響

    2.2.1 應(yīng)力-應(yīng)變關(guān)系

    不同凍融次數(shù)、不同圍壓下的應(yīng)力-應(yīng)變曲線如圖4 所示,其中F-T 表示凍融次數(shù)。從圖4 可以發(fā)現(xiàn),圍壓100 kPa 時應(yīng)力-應(yīng)變曲線呈現(xiàn)出應(yīng)變軟化型,經(jīng)過多次凍融循環(huán),有向應(yīng)變硬化型轉(zhuǎn)變的趨勢,剪切破壞時試樣產(chǎn)生破裂面;圍壓200、300 kPa 時表現(xiàn)出應(yīng)變硬化型,破壞形式也由脆性轉(zhuǎn)變?yōu)樗苄?,試樣發(fā)生剪脹破壞。通過不同凍融次數(shù)的應(yīng)力-應(yīng)變曲線可以看出,前期凍融循環(huán)作用對土體強度的影響十分敏感,尤其是第1 次和第3 次凍融結(jié)束后,不同圍壓下的應(yīng)力-應(yīng)變曲線均明顯低于其他凍融次數(shù),隨著凍融次數(shù)增加,強度有所增大并趨于穩(wěn)定。

    圖4 不同凍融次數(shù)、不同圍壓的應(yīng)力-應(yīng)變曲線Fig.4 Stress-strain curves of different freeze-thaw times and different confining pressures

    2.2.2 破壞強度

    凍融次數(shù)與破壞強度的關(guān)系如圖5 所示,可以看出,經(jīng)過3 次的凍融循環(huán)后土體結(jié)構(gòu)產(chǎn)生損傷,破壞強度急劇降低。根據(jù)圖3 試樣的凍縮變化率發(fā)現(xiàn),凍融1~3 次的高度凍縮率曲線的斜率是3~5 次的1.97 倍,凍融作用的影響明顯減弱。凍縮速率的快速減小代表內(nèi)部結(jié)構(gòu)將不會再進行劇烈調(diào)整,隨著土體的凍縮,固結(jié)度提高,結(jié)構(gòu)強度逐漸恢復;另外,第5 次凍融后的試樣抗剪強度發(fā)生大幅回升,說明第4 次凍融后土體抗剪強度應(yīng)增加而非降低。因此,凍融循環(huán)的劣化作用發(fā)生在前3 次。雖然多次凍融后強度趨于穩(wěn)定,但仍低于未凍融時的強度,說明凍融作用對土體造成的損傷具有不可恢復性。

    圖5 不同凍融次數(shù)下破壞強度變化Fig.5 Change of failure strength under different freeze-thaw times

    2.2.3 抗剪強度指標

    內(nèi)摩擦角和黏聚力在不同凍融次數(shù)下的變化規(guī)律如圖6 所示。內(nèi)摩擦角經(jīng)歷1 次凍融之后便達到平衡;黏聚力與強度變化相似,先減小后增加并達到穩(wěn)定。經(jīng)分析,初次凍融后,由于水分的相變和遷移,體積膨脹擠壓鄰近土顆粒,破壞土顆粒間的膠結(jié),降低咬合摩擦力,引起內(nèi)摩擦角和黏聚力均減小[22]。試樣的凍縮作用增大固結(jié)度,顆粒團聚體解體,小顆粒增多,咬合程度提高,所以內(nèi)摩擦角受后續(xù)的凍融影響較小[3]。土體結(jié)構(gòu)凍融過程中的水分運移、溶蝕、沖刷黏性礦物和產(chǎn)生的冷生結(jié)構(gòu)等削弱了顆粒間的膠結(jié)力,導致黏聚力繼續(xù)降低,隨著凍融次數(shù)增加,水分運移通道趨于穩(wěn)定,新的土體結(jié)構(gòu)逐漸形成,因此黏聚力也不再發(fā)生較大變化。

    圖6 不同凍融次數(shù)下抗剪強度參數(shù)變化Fig.6 The change of shear strength parameters under different freeze-thaw times

    3 理論分析與離散元數(shù)值模擬

    依據(jù)室內(nèi)試驗結(jié)果,本文將通過離散元數(shù)值模擬研究前3 次凍融循環(huán)作用引起強度損傷的細觀影響機理和作用過程。

    3.1 離散元單元間接觸模型

    顆粒離散元法中,巖土體是由一系列遵循牛頓運動定律的離散單元堆積而成,單元間依靠可斷裂的彈簧表示法向力(Fn)、切向力(Fs)和摩擦力(Ff)等效真實世界中巖石或土壤中礦物顆粒之間的膠結(jié)作用。接觸模型參數(shù)中Kn為法向剛度,Ks為切向剛度,Xb為斷裂位移,Xn為法向相對位移,Xs為切向相對位移,F(xiàn)s0和μp分別為切向初始抗剪力和摩擦因數(shù),如圖7 所示。

    圖7 離散元顆粒接觸關(guān)系Fig.7 Contact of discrete element

    單元間法向力與法向剛度的關(guān)系可用下式表示:

    相鄰單元連接時相互之間存在彈簧的拉壓作用,若位移Xn超過斷裂位移Xb時彈簧斷裂,單元間將不會產(chǎn)生拉力,但是發(fā)生擠壓時仍然可以產(chǎn)生壓力作用。

    單元間切向力與切向剛度的關(guān)系為:

    根據(jù)摩爾庫倫準則,單元間連接完整時最大剪切力Fsmax表示為:

    當剪切力超過Fsmax時,彈簧發(fā)生斷裂,初始抗剪力Fs0消失,發(fā)生滑動,產(chǎn)生的摩擦力Ff為-μpFn,與真實顆粒相比,膠結(jié)斷裂意味著法向和切向之間的連接均會斷開。

    由于在物理世界中巖土體的顆粒并不是完全彈性體,內(nèi)部存在水分、節(jié)理或微裂隙等因素阻礙顆粒的運動,動能會逐漸轉(zhuǎn)化為熱能而衰減。在數(shù)值模擬中,單元上的黏滯阻力(Fv)表征了真實世界中各種因素對動能的消耗,MatDEM 通過引入該力避免了動能在系統(tǒng)中的積累[23]。黏滯阻力由下式給定:

    式中:ψ為阻尼系數(shù),與模型的體積、顆粒直徑、顆粒質(zhì)量和顆粒剛度系數(shù)有關(guān);μ為顆粒的速度。

    3.2 離散元多場耦合理論

    在離散元數(shù)值模擬中,模型中的顆粒代表由一定體積黏土顆粒和水組成的集合體,其中水是虛設(shè)值,并不真實存在。顆粒除了具有一定的力學性質(zhì)外,還可以通過虛設(shè)一個溫度值,賦給顆粒后讓其具有溫度屬性。在本研究中,真實土體的條件是含水率為17.6%,溫度為20℃,將這2 種條件賦到顆粒上便可以通過建立半徑與溫度、含水率和應(yīng)力之間的關(guān)系,實現(xiàn)凍融循環(huán)中溫度場、水分場和應(yīng)力場的耦合模擬。假定顆粒的含水率表示土顆粒及其周圍區(qū)域的平均含水率,溫度和水分通過顆粒間的接觸進行傳遞。

    由于水凍結(jié)時會釋放潛熱出現(xiàn)過冷現(xiàn)象[24],模擬時假設(shè)將0℃至過冷溫度的區(qū)域設(shè)為凍結(jié)緣,且水分僅在該區(qū)域內(nèi)運移[25]。

    3.2.1 溫度場和水分場控制方程

    在數(shù)值模擬中,首先進行溫度傳遞,單元顆粒經(jīng)過一個時間步的溫度變化量?T,通過利用溫度場控制方程[26]轉(zhuǎn)化為下式進行計算:

    式中:C為土的體積熱容,kJ/(m3·K);λ為土與水的平均導熱系數(shù),W/(m·K);η為顆粒間溫度傳遞影響系數(shù),與時間步數(shù)和循環(huán)次數(shù)有關(guān);L為冰-水相變潛熱,一般取值為334.56 kJ/kg;T為單元顆粒初始溫度,℃;Tj為周圍接觸到的單元顆粒溫度,j=1,2,···,N;ρi為冰密度,取值0.917×103kg/m3;Vi為溫度低于0℃的顆粒單元體積;V為組成模型的單元顆??傮w積。

    隨著溫度的降低,顆粒間逐漸出現(xiàn)水分運移現(xiàn)象。每個時間步的水分運移量?ω參考水分場控制方程[27],轉(zhuǎn)換為如下公式進行計算,升溫時也按照此公式進行水分回遷模擬。

    式中:D為水分擴散系數(shù);ρw為水密度,取 值1×103kg/m3;ζ為顆粒溫度差與水分運移量的系數(shù),通過時間步數(shù)和循環(huán)次數(shù)等模擬訓練獲得;ωk為與周圍接觸到的單元顆粒的含水率,k=1,2,···,b;Tf為過冷溫度,取-1.5℃[28]。

    水分擴散系數(shù)D計算公式為:

    式中:a、m、l為水土特征曲線擬合參數(shù),引用基于室內(nèi)試驗數(shù)據(jù)擬合的V-G 模型和凍土物理學[29];θi為冰的體積含量;S為相對飽和度;θw為液態(tài)水的體積含量;ws為飽和含水率;wc為殘余含水率;Kb為飽和滲透系數(shù)。

    3.2.2 單元顆粒半徑變化及對強度的影響

    在凍融過程中,由于溫度的變化,土骨架凍縮、回彈和水分的運移及相變均會對土體體積產(chǎn)生影響,為了在數(shù)值模擬中體現(xiàn)這一過程,本文將上述變化細化到單元顆粒的變形上,把單元顆粒的變形原因視為土顆粒和水分變化兩部分。假設(shè)顆粒的半徑按照土顆粒和水分之間的質(zhì)量比值進行劃分,即水分所占半徑的大小為含水率的值ω,土顆粒占半徑的(1-ω)。對于土顆粒,計算時采用土顆粒所占半徑的比例乘以土顆粒收縮系數(shù),模擬土顆粒受溫度影響引起的顆粒半徑變化。對于水分,考慮了水分運移引起的含水率增加或減少以及水分凍結(jié)和融化導致的顆粒半徑脹縮變化。為了展示變化過程,給出了顆粒在含水率增加時的凍融過程圖如圖8 所示。

    圖8 含水率增加時顆粒變化Fig.8 Particle changes when the moisture content increases

    因此,當顆粒溫度低于過冷溫度或高于0℃時,顆粒半徑R的變化采用如下公式計算:

    式中:R0為初始單元顆粒半徑;R1為融化時單元顆粒半徑;ω0為初始含水率;ωd、ωr分別為凍結(jié)完成和融化完成后的含水率;γ為土顆粒收縮相關(guān)系數(shù),根據(jù)試樣凍縮量取值。溫度低于過冷溫度視為凍結(jié)完成,高于0℃時視為融化完成,顆粒半徑不變。

    依據(jù)應(yīng)力場基本方程中應(yīng)力σ與彈性模量E、泊松比v和體積變形有關(guān)[30],模型中顆粒間的應(yīng)力場計算公式如下:

    式中:RB為凍結(jié)前顆粒半徑;RL為凍結(jié)或融化后的顆粒半徑。

    由于顆粒半徑R和泊松比v與顆粒間的法向剛度和切向剛度之間存在關(guān)系,從Liu Chun 等[31]給出的參數(shù)轉(zhuǎn)換公式可知:

    通過式(1)、式(2)可知法向力和切向力與法向剛度和切向剛度相關(guān),所以模型中顆粒的正應(yīng)力 σn和切應(yīng)力σs的計算公式為:

    式中:FnB、FsB與FnL、FsL分別為顆粒半徑為RB和RL時的法向力和切向力;xnB、xnL為顆粒的法向相對位移;xsB、xsL為切向相對位移。

    根據(jù)室內(nèi)試驗結(jié)果,試樣強度劣化發(fā)生在前3 次的凍融循環(huán),將彈性模量、泊松比與凍融次數(shù)結(jié)合,擬合公式如下:

    式中:B為凍融次數(shù)。

    在數(shù)值模擬單次迭代中,首先利用式(5)實現(xiàn)顆粒間的熱量傳遞;利用式(6)實現(xiàn)水分在溫度的影響下發(fā)生運移;然后篩選溫度低于凍結(jié)溫度或高于0℃的顆粒,利用式(9)或式(10)分別實現(xiàn)凍結(jié)和融化時溫度、水分運移及相變對顆粒半徑產(chǎn)生的影響;顆粒半徑的增大和縮小引起顆粒間的擠壓和移動,改變了顆粒間的力學性質(zhì),通過式(14)或式(15)計算當前狀態(tài)下的應(yīng)力;最后平衡模型。由于顆粒位置發(fā)生重排列,顆粒間的接觸狀態(tài)被改變,進而影響下一次迭代中溫度傳遞、水分運移和應(yīng)力變化。經(jīng)過以上步驟的循環(huán)迭代實現(xiàn)凍融條件下土體內(nèi)部的水熱力三場耦合。

    3.3 離散元模型的建立

    3.3.1 堆積模型的建立

    本模擬試驗基于南京大學自主研發(fā)的三維離散元模擬軟件MatDEM[32]。通過該軟件首先建立一個內(nèi)徑3.91 cm 的圓管,在圓管內(nèi)生成平均半徑0.1 cm 的隨機顆粒,然后對單元顆粒施加隨機速度進行壓密沉積,建立直徑0.039 1 m、高0.08 m 的模型如圖9 所示。由圓管、頂板及底板構(gòu)成的恒溫結(jié)構(gòu)如圖10 所示,溫度設(shè)定為20℃或-20℃,不參與水分運移和不對模型產(chǎn)生約束,僅作為溫度邊界。模型共由24 466 個活動單元組成,其中試樣模型由17 930 個單元顆粒組成。

    圖9 試樣模型Fig.9 Sample model

    圖10 恒溫結(jié)構(gòu)Fig.10 Thermostatic structure

    3.3.2 模型的材料設(shè)置

    采用離散元法建模時,為了獲得與真實土體宏觀力學性質(zhì)相似的堆積模型,在MatDEM 軟件中的材料自動訓練模塊,利用室內(nèi)試驗獲取的土體力學參數(shù)如:彈性模量(E)、泊松比(v)、抗壓強度(Cu)、拉伸強度(Tu)和摩擦因數(shù)(μp),進行連續(xù)5 次的訓練優(yōu)化,獲得建立模型所需的宏觀力學參數(shù),然后根據(jù)Liu Chun 等[31]給出的宏微觀參數(shù)轉(zhuǎn)換公式,計算反映模型顆粒間性質(zhì)的微觀力學參數(shù),如:法向剛度(Kn)、剪切剛度(Ks)、初始抗剪力(Fs0)、斷裂位移(Xb)和摩擦因數(shù)(μp)。完成材料訓練后,將獲得的力學參數(shù)賦給堆積的試樣模型,進行單軸無側(cè)限壓縮模擬試驗,并與室內(nèi)無側(cè)限抗壓試驗結(jié)果對比,如圖11 所示。室內(nèi)試驗所獲的土體力學參數(shù)和通過對模型材料自動訓練獲得的宏微觀力學參數(shù)見表1。

    表1 土體室內(nèi)試驗力學參數(shù)和模型材料的宏微觀力學參數(shù)Table 1 Macroscopic parameters of soil test maro-micro parameters of model meterials

    圖11 室內(nèi)與數(shù)值模擬的無側(cè)限抗壓強度Fig.11 Indoor and numerical simulation of unconfined compressive strength

    數(shù)值模擬中溫度控制方程和水分傳導方程的相關(guān)參數(shù)見表2,其中η和 ζ經(jīng)過多次訓練分別取0.1 和0.05,水和冰的導熱系數(shù)λw和λi均取自規(guī)范值[29],式(7)、式(8)中ws、wc、Kb及式(5)中C、λ根據(jù)已有的室內(nèi)土工試驗成果進行取值。

    表2 分析計算的相關(guān)參數(shù)Table 2 Macroscopic and microscopic mechanical parameters of soil

    4 數(shù)值模擬結(jié)果分析

    4.1 溫度場演變規(guī)律

    試樣的凍融過程如圖12 所示。根據(jù)模型中心區(qū)域的3 個溫度節(jié)點:0℃、過冷溫度和-20℃,將凍結(jié)過程劃分為溫度快速下降、緩慢相變過程、繼續(xù)降溫、溫度穩(wěn)定4 個階段。第一階段,如圖12a-圖12e 所示,在50 min 時由于試樣表層與恒溫結(jié)構(gòu)之間溫差大,顆??焖賰鼋Y(jié),形成明顯凍結(jié)區(qū);凍結(jié)區(qū)內(nèi)側(cè)分布著溫度位于0℃附近的一層較薄區(qū)域,形成凍結(jié)緣區(qū)[33];距離冷端較遠的顆粒溫度未明顯下降,形成未凍土區(qū)。至 170 min 時,因凍結(jié)面與未凍區(qū)溫差變小,所以凍結(jié)區(qū)擴大較慢,但是內(nèi)部溫度降低較快,模型中心區(qū)域溫度基本全部降至0℃左右,結(jié)合圖12c 和圖12d 可以發(fā)現(xiàn),凍結(jié)鋒面推進較多,凍結(jié)緣變得越來越厚。凍結(jié)進入第二階段,如圖12f 和圖12g 所示,此時凍結(jié)緣內(nèi)冰水共存,水相變時釋放潛熱,減緩降溫速度。直到250 min 后模型核心溫度全部低于過冷溫度,可視為模型已完全凍結(jié)。第三階段開始,此時模型中溫度的傳遞速度有所提高,顆粒溫度繼續(xù)降低,到400 min后(圖12h),降溫過程結(jié)束,進入第四階段,顆粒與恒溫結(jié)構(gòu)沒有溫差,達到溫度穩(wěn)定狀態(tài)。

    圖12 凍融過程中的溫度傳遞Fig.12 Temperature transfer during freeze-thaw

    融化過程與凍結(jié)過程類似,同樣經(jīng)過4 個階段。將圖12i-圖12k 所處時間點的溫度變化與凍結(jié)過程相比,區(qū)別在于融化過程中0℃以下的變溫速率大于0℃以上的變溫速率,直至全部融化。第2、第3 次凍融循環(huán)中的溫度傳遞過程與第1 次基本一致。

    4.2 水分場演變規(guī)律

    依據(jù)4.1 節(jié)中對凍融過程的階段劃分,可以看出水分運移范圍跟隨凍結(jié)區(qū)向模型的內(nèi)部逐漸轉(zhuǎn)移。在凍結(jié)過程的第一階段中,如圖13a-圖13e 所示,由于冷端溫度與未凍區(qū)之間溫差大,土中水分成冰較快,水分運移量很少。第二階段中,如圖13f-圖13g 所示,凍結(jié)區(qū)發(fā)展緩慢,但凍結(jié)緣增厚,土中水分有較長的時間進行遷移,所以圖13e-圖13g 中含水率升高的顆粒明顯增多。第三階段中,模型已完全凍結(jié),含水率基本不變,如圖13h 所示。以上模擬水分運移過程說明,溫度梯度的存在導致引起水分運移的土水勢出現(xiàn),水分從高溫向低溫處遷移[34];溫差越大,土體凍結(jié)越快,水分以原位凍結(jié)為主;溫度梯度較小時,凍結(jié)較慢,水分運移增加。凍結(jié)區(qū)邊緣顆粒的體積會通過周圍顆粒的水分補給逐漸增大,擠壓周圍顆粒,說明水分運移會引起水分重分布,局部含水率升高,水分凍結(jié)成冰后導致土體發(fā)生不均勻變形,改變顆粒排列方式,破壞內(nèi)部結(jié)構(gòu)。

    圖13 凍融過程中水分運移Fig.13 Moisture migration during freeze-thaw

    融化過程中,在溫度梯度和土水勢的作用下,土體水分發(fā)生回遷現(xiàn)象,高含水率顆粒中的水分向周圍低含水率的顆粒轉(zhuǎn)移,如圖13i-圖13l 所示。升溫結(jié)束后,發(fā)現(xiàn)水分并未完全回遷,導致含水率分布不均,說明顆粒凍結(jié)后位置被動調(diào)整,與周圍顆粒的連接狀態(tài)發(fā)生改變,所以不能按照水分運移的路徑回遷。在第2、第3 次凍融循環(huán)中,水分的遷移會繼續(xù)引起土體結(jié)構(gòu)的變化,含水率升高或降低的顆粒逐漸增加,進而對溫度傳遞、水分運移的路徑和結(jié)構(gòu)強度產(chǎn)生影響。

    4.3 單元顆粒半徑變化及宏觀變形

    圖14 為17 930 個模型顆粒半徑在首次凍結(jié)后和融化后的大小分布圖,可以發(fā)現(xiàn)顆粒半徑主要集中在0.000 75~0.001 10 m,凍融后的平均半徑低于0.001 m,發(fā)生細觀角度上土顆粒和水分集合體的凍縮現(xiàn)象。受水分運移的影響,凍融前后相同半徑區(qū)間的單元顆粒數(shù)均有差別,融化后的分布曲線偏低,顆粒大小更為分散。隨著凍融次數(shù)增加,含水率升高和降低的顆粒增多,最大半徑和最小半徑的顆粒也隨之增長。顆粒半徑的變化引起周圍顆粒發(fā)生位移,原本的接觸狀態(tài)隨之改變,形成位移場,如圖15 中顆粒在X方向的位移。通過圖15a、圖15b 對比發(fā)現(xiàn),前后位移場存在一定的差異,說明模型在凍融作用下內(nèi)部結(jié)構(gòu)受顆粒的脹縮影響已經(jīng)發(fā)生變化。

    圖14 凍結(jié)、融化后顆粒半徑分布Fig.14 Particle radius distribution after freezing and thawing

    圖15 凍結(jié)、融化后X 方向位移場Fig.15 X-directional displacement field after freezing and thawing

    通過對模型凍結(jié)和融化后的直徑和高度變化數(shù)據(jù)的提取,凍結(jié)后含水率升高的顆粒半徑明顯增加,擠壓周圍顆粒,引起模型凍脹;融化后,顆粒半徑縮小,在顆粒間膠結(jié)力、黏聚力等的作用下向內(nèi)側(cè)移動,變形逐漸恢復(圖16)。整體上可以看出模型在凍融過程中體積先凍脹而后融沉的現(xiàn)象。由于顆粒代表一定體積的土水集合體,顆粒直徑占高度的1/80,直徑的1/39,所以宏觀變形上較難達到一致,但是模擬的最高誤差未超過0.172%,可以有效說明凍融作用對土體變形的影響。

    圖16 凍結(jié)、融化后高度和直徑變化率Fig.16 Rate of change in height and diameter after freezing and thawing

    4.4 溫度和水分綜合作用對應(yīng)力和強度的影響

    凍融作用引起的內(nèi)部應(yīng)力變化如圖17 所示。凍結(jié)前顆粒間的應(yīng)力受重力影響由上至下逐漸增加,但是變化并不明顯。經(jīng)過凍融作用,在溫度、水分運移及相變的影響下引起部分顆粒的脹縮讓模型內(nèi)產(chǎn)生局部的變形,引起顆粒重新排列,導致顆粒斷開連接或與其他顆粒建立新連接,改變了顆粒間的受力狀態(tài),破壞了土體的內(nèi)部結(jié)構(gòu),最終造成應(yīng)力明顯降低。

    圖17 凍結(jié)前、融化后X 方向應(yīng)力場Fig.17 X-direction stress field before freezing and after melting

    將經(jīng)過0、1 和3 次凍融過程的模型進行圍壓為100 kPa 的單軸壓縮模擬試驗,與室內(nèi)相同條件的試驗結(jié)果進行對比,結(jié)果如圖18 所示,模型破壞強度大小與室內(nèi)試驗結(jié)果基本相同。在凍融循環(huán)模擬中,顆粒的脹縮、冰-水相變、水分聚集與回遷現(xiàn)象的發(fā)生,改變了土體原有的穩(wěn)定結(jié)構(gòu)。隨著凍融次數(shù)的增加,顆粒之間的孔隙、位置、連接狀態(tài)和應(yīng)力大小不斷變化。其次,凍結(jié)時產(chǎn)生的水分運移令部分顆粒持續(xù)獲得水分而明顯增大,與土體中局部孔隙得到水的補充后,凍結(jié)成冰發(fā)生膨脹現(xiàn)象一致,改變了顆粒間的相對位置及接觸狀態(tài),最終造成難以恢復的結(jié)構(gòu)損傷,導致土體強度降低。

    圖18 不同凍融次數(shù)的強度對比Fig.18 Strength comparison of different freeze-thaw crcles

    利用MatDEM 軟件中的材料訓練模塊和能對單元顆粒添加虛設(shè)值的特點,獲得與真實土體溫度和含水率相同、力學性質(zhì)相似的模型,并通過溫度與含水率變化影響顆粒半徑,引起的模型強度和變形與室內(nèi)試驗結(jié)果基本一致,說明數(shù)值模擬的結(jié)果是可靠的,且具有實際參考價值。

    5 結(jié)論

    a.土體在凍融循環(huán)作用下發(fā)生凍縮,體積減??;抗剪強度及其參數(shù)的劣化主要受前3 次凍融作用的影響,第3 次劣化程度最大,隨著凍融次數(shù)的增加,抗剪強度及其參數(shù)回升并趨于穩(wěn)定。

    b.建立的水、熱、力和半徑變化的計算方程,適用于非飽和黏性土;可較好地模擬溫度、水分和應(yīng)力的動態(tài)變化規(guī)律,獲得的試樣變形和破壞強度與室內(nèi)試驗結(jié)果吻合度也較高;可用于描述凍結(jié)鋒面的移動軌跡和含水率的變化過程。

    c.溫度傳遞過程中可劃分為溫度快速下降、緩慢相變過程、繼續(xù)降溫、溫度穩(wěn)定4 個階段;水分運移主要發(fā)生在溫度傳遞的前2 個階段,緩慢相變過程中的水分運移量居多。

    d.凍融作用導致顆粒整體縮小且分散性增大;土體主要由于土顆粒的脹縮、冰-水相變、水分運移和冷生結(jié)構(gòu)等共同作用引起位移場和應(yīng)力場改變,產(chǎn)生不可逆的結(jié)構(gòu)損傷而造成抗剪強度劣化。

    猜你喜歡
    模型
    一半模型
    一種去中心化的域名服務(wù)本地化模型
    適用于BDS-3 PPP的隨機模型
    提煉模型 突破難點
    函數(shù)模型及應(yīng)用
    p150Glued在帕金森病模型中的表達及分布
    函數(shù)模型及應(yīng)用
    重要模型『一線三等角』
    重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
    3D打印中的模型分割與打包
    大又大粗又爽又黄少妇毛片口| 99久久九九国产精品国产免费| 最近中文字幕2019免费版| 亚洲av国产av综合av卡| 欧美少妇被猛烈插入视频| 免费观看的影片在线观看| 免费看日本二区| 欧美潮喷喷水| 永久免费av网站大全| 国产成人免费观看mmmm| 免费观看性生交大片5| 久久久久国产精品人妻一区二区| 免费看不卡的av| 午夜福利在线观看免费完整高清在| 搡女人真爽免费视频火全软件| 综合色丁香网| 一级av片app| 日韩在线高清观看一区二区三区| 2021天堂中文幕一二区在线观| 色哟哟·www| 99久久精品国产国产毛片| 国产在线一区二区三区精| 成人高潮视频无遮挡免费网站| 汤姆久久久久久久影院中文字幕| 深爱激情五月婷婷| 女的被弄到高潮叫床怎么办| 少妇高潮的动态图| 99热网站在线观看| 毛片女人毛片| 日韩一区二区视频免费看| 国产免费福利视频在线观看| 国国产精品蜜臀av免费| 国产精品国产三级国产专区5o| 人妻 亚洲 视频| 国产女主播在线喷水免费视频网站| 最近最新中文字幕免费大全7| 午夜福利视频精品| 精品人妻一区二区三区麻豆| av.在线天堂| www.色视频.com| av卡一久久| 亚洲人成网站在线播| 国产一区二区三区综合在线观看 | 国产国拍精品亚洲av在线观看| 午夜福利视频1000在线观看| 免费黄频网站在线观看国产| 国产免费一级a男人的天堂| 日韩三级伦理在线观看| 国产亚洲精品久久久com| 观看免费一级毛片| 国产熟女欧美一区二区| 最近最新中文字幕免费大全7| 欧美xxxx黑人xx丫x性爽| 免费黄网站久久成人精品| 黄色配什么色好看| 亚洲精品成人久久久久久| 黄色怎么调成土黄色| 免费看av在线观看网站| 大香蕉久久网| 又大又黄又爽视频免费| 精品久久久精品久久久| 成人美女网站在线观看视频| 麻豆久久精品国产亚洲av| 亚洲国产欧美在线一区| 久久人人爽人人片av| 久久这里有精品视频免费| 交换朋友夫妻互换小说| 欧美一区二区亚洲| 蜜桃久久精品国产亚洲av| 99久久九九国产精品国产免费| 成人国产麻豆网| 久久韩国三级中文字幕| 观看免费一级毛片| 制服丝袜香蕉在线| 在线观看国产h片| 国产精品一二三区在线看| 久久亚洲国产成人精品v| 久久久午夜欧美精品| 搡老乐熟女国产| 韩国av在线不卡| 国产精品人妻久久久久久| freevideosex欧美| 久久99热这里只频精品6学生| 国产精品伦人一区二区| 老司机影院毛片| 亚洲av男天堂| 六月丁香七月| 黄色欧美视频在线观看| 欧美亚洲 丝袜 人妻 在线| 中国国产av一级| 热re99久久精品国产66热6| 欧美人与善性xxx| 哪个播放器可以免费观看大片| 亚洲精品久久午夜乱码| 亚洲国产精品国产精品| 22中文网久久字幕| 久久精品久久精品一区二区三区| 又爽又黄a免费视频| av福利片在线观看| 国产av不卡久久| 国产黄色视频一区二区在线观看| 国产成人freesex在线| eeuss影院久久| 亚洲国产精品成人综合色| av线在线观看网站| 午夜福利网站1000一区二区三区| 一级二级三级毛片免费看| 免费看av在线观看网站| 人体艺术视频欧美日本| 色哟哟·www| 久久久久性生活片| 人人妻人人澡人人爽人人夜夜| 少妇人妻久久综合中文| 亚洲成人久久爱视频| 国产又色又爽无遮挡免| 久久精品国产自在天天线| 欧美日韩国产mv在线观看视频 | 亚洲第一区二区三区不卡| 自拍欧美九色日韩亚洲蝌蚪91 | 精品久久久精品久久久| 人人妻人人看人人澡| 日韩欧美精品v在线| 国产av国产精品国产| 高清在线视频一区二区三区| 国产乱人视频| 日日摸夜夜添夜夜爱| 中文在线观看免费www的网站| 亚洲欧美精品自产自拍| 18禁在线播放成人免费| 最近最新中文字幕大全电影3| 在线看a的网站| 亚洲人成网站高清观看| 亚洲欧美精品自产自拍| 校园人妻丝袜中文字幕| 国产成人精品一,二区| 国产爱豆传媒在线观看| 亚洲,一卡二卡三卡| 精品人妻一区二区三区麻豆| 国产精品国产三级专区第一集| 亚洲精品久久午夜乱码| 成人综合一区亚洲| 97人妻精品一区二区三区麻豆| 一个人观看的视频www高清免费观看| 伊人久久国产一区二区| 日产精品乱码卡一卡2卡三| 精品少妇久久久久久888优播| 久久精品国产a三级三级三级| 肉色欧美久久久久久久蜜桃 | 亚洲国产日韩一区二区| 一级片'在线观看视频| 一级黄片播放器| 国产成人freesex在线| 亚洲av.av天堂| 水蜜桃什么品种好| 高清视频免费观看一区二区| 亚洲自偷自拍三级| 国产v大片淫在线免费观看| .国产精品久久| 国产亚洲91精品色在线| 亚洲国产成人一精品久久久| 两个人的视频大全免费| 一二三四中文在线观看免费高清| 中文字幕人妻熟人妻熟丝袜美| 日本黄大片高清| 国产亚洲最大av| 欧美日韩亚洲高清精品| 中文字幕av成人在线电影| 亚洲av福利一区| 国模一区二区三区四区视频| 六月丁香七月| 国产老妇伦熟女老妇高清| 亚洲国产色片| 中文字幕亚洲精品专区| tube8黄色片| 波野结衣二区三区在线| 少妇人妻 视频| 国产综合精华液| 国产 精品1| 老师上课跳d突然被开到最大视频| 久久女婷五月综合色啪小说 | 色视频在线一区二区三区| 午夜激情福利司机影院| 久久久欧美国产精品| 国产午夜精品久久久久久一区二区三区| 国产黄片视频在线免费观看| 干丝袜人妻中文字幕| 亚洲色图av天堂| 亚洲精品乱码久久久久久按摩| 精品一区在线观看国产| 国产v大片淫在线免费观看| 网址你懂的国产日韩在线| 男人舔奶头视频| 久久久国产一区二区| 禁无遮挡网站| 美女被艹到高潮喷水动态| 久久99蜜桃精品久久| 久久影院123| 日韩一区二区三区影片| 国产精品一区二区三区四区免费观看| 亚洲国产av新网站| 嫩草影院入口| 国产在线一区二区三区精| 国产免费一级a男人的天堂| 久久影院123| 中文字幕亚洲精品专区| 精品国产露脸久久av麻豆| 亚洲国产欧美人成| 免费观看的影片在线观看| 男女那种视频在线观看| 男女啪啪激烈高潮av片| 一边亲一边摸免费视频| 啦啦啦啦在线视频资源| 成人免费观看视频高清| xxx大片免费视频| 一级黄片播放器| 成人高潮视频无遮挡免费网站| 亚洲真实伦在线观看| 欧美精品国产亚洲| 人妻系列 视频| 亚洲色图av天堂| 18禁裸乳无遮挡动漫免费视频 | 五月开心婷婷网| 精品国产三级普通话版| 美女cb高潮喷水在线观看| 美女被艹到高潮喷水动态| 国产 一区精品| 日本一二三区视频观看| 久久久a久久爽久久v久久| 久久久久久久久久久免费av| 亚洲伊人久久精品综合| 18禁裸乳无遮挡免费网站照片| 蜜桃亚洲精品一区二区三区| 另类亚洲欧美激情| 免费观看在线日韩| 中国国产av一级| 嘟嘟电影网在线观看| 在线天堂最新版资源| 男女国产视频网站| 如何舔出高潮| 1000部很黄的大片| 久久人人爽人人爽人人片va| 一区二区三区乱码不卡18| 最近2019中文字幕mv第一页| 一级av片app| 国产在线男女| 午夜福利高清视频| 久久精品综合一区二区三区| 欧美xxxx性猛交bbbb| 久热这里只有精品99| 婷婷色麻豆天堂久久| 久久久精品欧美日韩精品| 免费看日本二区| 少妇高潮的动态图| 亚洲av男天堂| 老司机影院成人| 美女被艹到高潮喷水动态| 亚洲av国产av综合av卡| 成人漫画全彩无遮挡| 观看免费一级毛片| 99久久精品热视频| 天堂俺去俺来也www色官网| 制服丝袜香蕉在线| 一区二区三区乱码不卡18| 搞女人的毛片| 久久ye,这里只有精品| 男女那种视频在线观看| 美女国产视频在线观看| 中文字幕制服av| 2022亚洲国产成人精品| 狂野欧美激情性bbbbbb| 亚洲av欧美aⅴ国产| 美女内射精品一级片tv| 嘟嘟电影网在线观看| 永久网站在线| 一级毛片我不卡| 亚洲在久久综合| 亚洲国产欧美人成| 久久精品熟女亚洲av麻豆精品| 国产精品av视频在线免费观看| 性色av一级| 91久久精品电影网| 久久久色成人| 永久网站在线| 99久久精品国产国产毛片| 白带黄色成豆腐渣| 久久6这里有精品| 中文资源天堂在线| 我要看日韩黄色一级片| 少妇的逼水好多| 日韩三级伦理在线观看| 丝袜喷水一区| 少妇丰满av| 国产老妇伦熟女老妇高清| 精品国产一区二区三区久久久樱花 | 亚洲激情五月婷婷啪啪| 国产淫片久久久久久久久| 美女脱内裤让男人舔精品视频| 久久6这里有精品| 亚洲无线观看免费| 18禁在线播放成人免费| 亚洲国产精品999| 日韩av免费高清视频| 国产大屁股一区二区在线视频| eeuss影院久久| 美女高潮的动态| 久久精品人妻少妇| 久久精品夜色国产| 日韩人妻高清精品专区| 亚洲欧美一区二区三区国产| 日韩人妻高清精品专区| 狂野欧美激情性bbbbbb| 男女边摸边吃奶| 男人爽女人下面视频在线观看| 亚洲成人av在线免费| 69av精品久久久久久| 国产成人一区二区在线| 亚洲av一区综合| 少妇的逼水好多| 国产男女内射视频| 久久精品熟女亚洲av麻豆精品| 丰满人妻一区二区三区视频av| av线在线观看网站| www.色视频.com| 国产爽快片一区二区三区| 欧美日韩视频精品一区| 男人狂女人下面高潮的视频| 91aial.com中文字幕在线观看| 久久人人爽人人片av| 在线观看av片永久免费下载| 亚洲成人久久爱视频| 麻豆久久精品国产亚洲av| 高清欧美精品videossex| 亚洲成人av在线免费| 毛片一级片免费看久久久久| 亚洲av二区三区四区| 国产黄片视频在线免费观看| 看免费成人av毛片| 卡戴珊不雅视频在线播放| 不卡视频在线观看欧美| 亚洲精品久久午夜乱码| 国产精品99久久久久久久久| 91在线精品国自产拍蜜月| 亚洲在久久综合| 国产综合精华液| av在线亚洲专区| 国产视频首页在线观看| 欧美日本视频| 青青草视频在线视频观看| 国产精品国产av在线观看| 欧美国产精品一级二级三级 | 伊人久久国产一区二区| 色吧在线观看| 男女那种视频在线观看| 国产成人精品婷婷| 80岁老熟妇乱子伦牲交| 听说在线观看完整版免费高清| 亚洲人成网站在线观看播放| 亚洲av电影在线观看一区二区三区 | 日日撸夜夜添| 久久久久久久亚洲中文字幕| xxx大片免费视频| 免费少妇av软件| 丰满人妻一区二区三区视频av| videos熟女内射| 国产精品久久久久久精品电影| 久久久久久久久久久丰满| 精品久久久久久电影网| 久久久a久久爽久久v久久| 亚洲一级一片aⅴ在线观看| 国产大屁股一区二区在线视频| 亚洲av男天堂| 天堂网av新在线| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 国产精品福利在线免费观看| 美女视频免费永久观看网站| 免费播放大片免费观看视频在线观看| av黄色大香蕉| 国产精品人妻久久久影院| 精品久久久噜噜| 精品午夜福利在线看| 国产一区有黄有色的免费视频| 亚洲精品色激情综合| 国产色爽女视频免费观看| 亚洲av免费在线观看| 亚洲在线观看片| 国产精品人妻久久久影院| 中文字幕久久专区| 国产亚洲最大av| 欧美日本视频| 亚洲欧美日韩卡通动漫| 欧美最新免费一区二区三区| 美女cb高潮喷水在线观看| 日本黄大片高清| 亚洲国产精品成人久久小说| 亚洲av国产av综合av卡| 国产精品伦人一区二区| 亚洲国产av新网站| 国产日韩欧美在线精品| 国产精品三级大全| 欧美老熟妇乱子伦牲交| 国产黄片视频在线免费观看| 免费av不卡在线播放| 国语对白做爰xxxⅹ性视频网站| 亚洲成人av在线免费| 婷婷色综合www| 国产免费一级a男人的天堂| 国产精品久久久久久精品电影小说 | 91aial.com中文字幕在线观看| 蜜桃亚洲精品一区二区三区| 精品久久久久久电影网| 国产免费又黄又爽又色| av一本久久久久| 色5月婷婷丁香| av国产精品久久久久影院| 国产毛片a区久久久久| 久久韩国三级中文字幕| 九九在线视频观看精品| 精品一区二区三区视频在线| 国产亚洲5aaaaa淫片| 亚洲成人一二三区av| 五月伊人婷婷丁香| av在线app专区| 亚洲国产欧美在线一区| 免费av观看视频| 午夜免费鲁丝| 亚洲成人av在线免费| 国内精品宾馆在线| 日韩一本色道免费dvd| 日韩一区二区三区影片| 在线观看一区二区三区| videos熟女内射| 一级片'在线观看视频| 国产精品国产av在线观看| 视频中文字幕在线观看| 国产高清不卡午夜福利| 国产精品嫩草影院av在线观看| 午夜福利高清视频| 国产精品伦人一区二区| 99热网站在线观看| 秋霞伦理黄片| 亚洲av一区综合| 看十八女毛片水多多多| 国产亚洲av嫩草精品影院| 交换朋友夫妻互换小说| 成人鲁丝片一二三区免费| 视频区图区小说| 亚洲人成网站高清观看| 亚洲最大成人av| 国产成人freesex在线| 好男人视频免费观看在线| 大话2 男鬼变身卡| 日韩一区二区三区影片| 日韩不卡一区二区三区视频在线| 天堂网av新在线| 国产精品爽爽va在线观看网站| 免费看光身美女| 丰满少妇做爰视频| 中文天堂在线官网| 男人狂女人下面高潮的视频| 一级毛片 在线播放| 精品视频人人做人人爽| 午夜免费鲁丝| 一级毛片黄色毛片免费观看视频| 最近的中文字幕免费完整| 成人国产麻豆网| 女的被弄到高潮叫床怎么办| 草草在线视频免费看| 久久久久久久亚洲中文字幕| 欧美性猛交╳xxx乱大交人| 五月开心婷婷网| 日本欧美国产在线视频| 国产精品成人在线| 97在线人人人人妻| 成人特级av手机在线观看| 欧美激情国产日韩精品一区| 精品久久久精品久久久| 亚洲精品乱码久久久久久按摩| 欧美 日韩 精品 国产| 少妇人妻久久综合中文| 午夜视频国产福利| 少妇猛男粗大的猛烈进出视频 | 性色av一级| 久久99热这里只有精品18| 黄色欧美视频在线观看| 午夜亚洲福利在线播放| 日本免费在线观看一区| 国产中年淑女户外野战色| 看非洲黑人一级黄片| 观看免费一级毛片| 青春草亚洲视频在线观看| 精品酒店卫生间| 一本一本综合久久| 少妇人妻一区二区三区视频| av在线老鸭窝| 免费大片18禁| 夜夜看夜夜爽夜夜摸| 免费在线观看成人毛片| 亚洲国产av新网站| 国产在线男女| 亚洲高清免费不卡视频| 日韩一本色道免费dvd| 亚洲国产成人一精品久久久| 午夜激情久久久久久久| 欧美+日韩+精品| 国语对白做爰xxxⅹ性视频网站| 亚洲精品亚洲一区二区| 高清毛片免费看| 亚洲国产av新网站| 国产乱人偷精品视频| 亚洲av免费在线观看| 日韩成人伦理影院| 免费观看性生交大片5| 亚洲最大成人中文| 91久久精品国产一区二区成人| 国产免费福利视频在线观看| 男人狂女人下面高潮的视频| 国产精品人妻久久久影院| 国产伦在线观看视频一区| 插阴视频在线观看视频| 亚洲精品,欧美精品| 欧美高清性xxxxhd video| 国产黄a三级三级三级人| 欧美人与善性xxx| 婷婷色综合www| 国产精品福利在线免费观看| 丝袜喷水一区| 综合色丁香网| 大话2 男鬼变身卡| 国产精品蜜桃在线观看| 午夜激情久久久久久久| 黄色一级大片看看| 黄色视频在线播放观看不卡| 欧美bdsm另类| 国产淫语在线视频| 国产免费视频播放在线视频| 欧美精品人与动牲交sv欧美| 成人亚洲精品av一区二区| 一边亲一边摸免费视频| 五月开心婷婷网| 黄色视频在线播放观看不卡| 一区二区av电影网| 另类亚洲欧美激情| 亚洲婷婷狠狠爱综合网| 亚洲精品自拍成人| 高清在线视频一区二区三区| 亚洲av电影在线观看一区二区三区 | 午夜亚洲福利在线播放| 国内精品美女久久久久久| 欧美成人精品欧美一级黄| 777米奇影视久久| 久久久久国产网址| 久久久久精品性色| 亚洲欧洲国产日韩| a级一级毛片免费在线观看| 六月丁香七月| 汤姆久久久久久久影院中文字幕| 在线观看人妻少妇| 亚洲精品久久午夜乱码| 亚洲av欧美aⅴ国产| 亚洲美女搞黄在线观看| 男女国产视频网站| 国产有黄有色有爽视频| 国产69精品久久久久777片| 亚洲av免费在线观看| 亚洲最大成人中文| 国产精品99久久99久久久不卡 | 日韩电影二区| 欧美成人a在线观看| 国产黄色免费在线视频| 丝袜喷水一区| 国产高清有码在线观看视频| 国产高潮美女av| 丝瓜视频免费看黄片| 搡老乐熟女国产| 成年女人看的毛片在线观看| 亚洲av不卡在线观看| 一级毛片电影观看| 熟女电影av网| 欧美成人精品欧美一级黄| 国国产精品蜜臀av免费| 美女脱内裤让男人舔精品视频| 欧美3d第一页| 国产黄a三级三级三级人| 免费av不卡在线播放| 韩国高清视频一区二区三区| 亚洲国产欧美在线一区| 久久久久久久久久久免费av| 国产视频内射| 大码成人一级视频| 成人免费观看视频高清| 99热这里只有是精品在线观看| 免费不卡的大黄色大毛片视频在线观看| 日本熟妇午夜| 久久久色成人| 一区二区av电影网| 国产伦在线观看视频一区| 国产淫片久久久久久久久| 久久99热这里只有精品18| 欧美+日韩+精品| 国产 精品1| 免费电影在线观看免费观看| 成人高潮视频无遮挡免费网站| 少妇 在线观看| 国产精品人妻久久久影院| 成年版毛片免费区| 少妇熟女欧美另类| 亚洲精品色激情综合| 白带黄色成豆腐渣| 性色av一级| 国产精品三级大全| 一级毛片黄色毛片免费观看视频| 欧美 日韩 精品 国产| 国产亚洲最大av| 日本黄大片高清| 久久精品国产鲁丝片午夜精品|