馬冬冬，汪鑫鵬，馬芹永1,，周志偉，楊 毅，袁 璞1,

（1.安徽理工大學深部煤礦采動響應與災害防控國家重點實驗室，安徽 淮南 232001；2.安徽理工大學土木建筑學院，安徽 淮南 232001；3.中國科學院西北生態(tài)環(huán)境資源研究院凍土工程國家重點實驗室，甘肅 蘭州 730000）

我國西部寒區(qū)以其獨特的地域特征和軍事作用在國家安全中處于非常重要的戰(zhàn)略地位[1]，在深厚凍土層以下，已建成大量以地下防護工程、深埋管道和地下戰(zhàn)略物資儲備庫為代表的國防工程[2-3]。在工程服役期間，需要考慮地面和鉆地武器對地下建筑物的爆破毀傷作用[4-6]，武器爆破后產生的應力波是結構損傷破壞的主要誘因[7]，研究凍土的動態(tài)本構關系是分析應力波在凍土內傳播規(guī)律的重要基礎和必要前提[8-9]，對凍土層以下防爆結構的優(yōu)化設計和安全穩(wěn)定性分析具有重要的參考價值。

目前，針對凍土動態(tài)力學性能和本構關系的研究多集中于單軸受力狀態(tài)[10-12]，但通過對凍土進行受力分析可知，在受到沖擊應力波作用前凍土已處于圍壓受力狀態(tài)，且隨著埋深的增加，凍土承受的限制作用逐漸增強，即所受的圍壓逐漸增大[13-15]，因此，需要研究圍壓作用下凍土的動態(tài)力學響應特征。馬芹永等[16]采用厚壁鋁質套筒限制凍土試樣的側向變形，以模擬被動圍壓受力狀態(tài)，利用分離式Hopkinson壓桿（split Hopkinson pressure bar, SHPB）試驗系統(tǒng)，對比了單軸和被動圍壓狀態(tài)下凍土的強度和變形破壞特征，發(fā)現(xiàn)單軸狀態(tài)下凍土的應力-應變曲線在達到峰值后逐漸下降，表現(xiàn)出應變軟化特征，而被動圍壓狀態(tài)下應力-應變曲線呈現(xiàn)明顯的應變硬化特征。Zhang 等[17]研究發(fā)現(xiàn)被動圍壓的存在能夠顯著提高凍土的動態(tài)峰值應力，并建立了能夠考慮應變率效應和溫度效應的凍土唯象本構模型，該模型能夠模擬圍壓對凍土應力-應變曲線峰前上升段的影響，但對曲線下降段的擬合效果有待進一步提高。馬冬冬等[18-19]利用能夠施加主動圍壓的液壓系統(tǒng)，開展了不同主動圍壓等級下凍結黏土和凍結砂土的SHPB 試驗，發(fā)現(xiàn)凍土在單軸狀態(tài)下呈脆性破壞特征，而主動圍壓狀態(tài)下凍土試樣呈塑性破壞。主動圍壓狀態(tài)下，凍結黏土和凍結砂土動態(tài)應力-應變曲線表現(xiàn)出不同的變形特征。圖1 為相同溫度和主動圍壓等級下凍結黏土和凍結砂土的動態(tài)應力-應變曲線，可以看出，凍結黏土的動態(tài)應力-應變曲線可分為彈性、塑性和破壞階段，且3 階段的過渡較平滑。相較于凍結黏土，凍結砂土達到塑性變形狀態(tài)所需的變形較小，即其彈性變形量遠小于凍結黏土，導致曲線中塑性階段的占比明顯增大，且彈性階段和塑性階段之間的拐點非常明顯。圍壓的存在會導致凍結黏土和凍結砂土的動態(tài)抗壓強度明顯增大，例如：當應變率為210～220 s-1、溫度為-15 ℃時，在1.5 MPa 主動圍壓下，凍結黏土和凍結砂土的動態(tài)抗壓強度相較于無圍壓狀態(tài)分別提高了61%和81%。此外，隨著圍壓等級的提高，凍結黏土的動態(tài)抗壓強度呈非線性增大，而凍結砂土的動態(tài)抗壓強度基本呈線性增大。相同圍壓條件下，凍結黏土和凍結砂土的峰值應力、峰值應變和極限應變均隨應變率的提高而增大?；诮浀涞闹?王-唐（ZWT）模型，通過在模型中引入被視為Weibull 分布的損傷變量和能夠考慮圍壓影響的應力增強因子，建立了主動圍壓狀態(tài)下凍結黏土的動態(tài)本構模型[20]。通過對比發(fā)現(xiàn)，該動態(tài)本構模型能夠較好地預測凍結黏土的動態(tài)應力-應變曲線特征和圍壓對其動態(tài)抗壓強度的增強效應。

圖1 主動圍壓作用下凍結黏土和凍結砂土的動態(tài)應力-應變曲線[18-19]Fig.1 Dynamic stress-strain curves of frozen clay and frozen sandy soil under active confining pressure[18-19]

綜合以上分析可知，作為一種典型的非線性黏彈性本構模型，ZWT 模型能夠較好地模擬凍結黏土的動態(tài)強度和變形特性，但無法較好地擬合主動圍壓作用下凍結砂土的動態(tài)應力-應變曲線特征。本文中，通過對經典的ZWT 模型進行改進，在模型中引入塑性體元件，并結合損傷力學理論，建立能夠考慮圍壓效應的凍結砂土動態(tài)損傷本構模型，通過將模型和試驗獲得的不同負溫和主動圍壓等級下凍結砂土的應力-應變曲線進行對比，對模型的預測效果進行驗證和評價。

1 考慮塑性變形的ZWT 模型

經典的ZWT 模型[21]由1 個非線性體和2 個Maxwell 體并聯(lián)而成，并通過控制松弛時間，反映不同應變率范圍內Maxwell 體的響應特征：

式中：E0、α 和β 均為與材料性質有關的彈性常數(shù)，E1和E2分別為低頻和高頻Maxwell 體的彈性常數(shù)，θ1和θ2分別為低頻和高頻Maxwell 體的松弛時間。

在沖擊荷載作用下，凍土以102～103s-1的應變率從損傷狀態(tài)迅速進入破壞狀態(tài)，特征時間為0.1 s 的低頻Maxwell 體無充足的松弛時間，將失去其作用并轉化為彈性體[22-23]，模型表達式可轉化為：

由圖1 可知，當應力超過一定值后，凍結砂土的動態(tài)應力-應變曲線將從線彈性階段迅速進入塑性階段，說明彈、塑性拐點（即屈服點）非常明顯，且塑性段會持續(xù)較長時間。因此，在模型的非線性體上串聯(lián)一個塑性體，即可得到能夠考慮塑性變形的ZWT 模型，見圖2，σs為塑性體的屈服強度。

圖2 考慮塑性變形的ZWT 模型Fig.2 The ZWT model established in consideration of plastic deformation

在非線性體和塑性體串聯(lián)后，模型表達式為：

沖擊荷載作用下，凍土試樣內部的原生裂紋和孔洞不斷擴展，并伴隨著新裂紋的產生，其損傷是一個逐漸累計的過程。凍土的動態(tài)損傷演化過程服從Weibull 分布規(guī)律[24]，即損傷變量表達式為：

式中：μ和η 為凍結砂土的損傷參數(shù)。

由式(3)可以看出，模型以屈服強度為分界點，可分為彈塑性屈服點之前和之后2 個部分。在屈服點前的彈性階段，凍結砂土試樣內部主要為彈性變形，而屈服點后的塑性階段，塑性變形導致試樣內部損傷逐步增大。因此，將損傷變量D代入式(3)中屈服點之后的表達式，即可得到綜合考慮損傷演化和塑性變形的ZWT 模型：

將動態(tài)應力-應變曲線中彈塑性拐點對應的縱坐標視為凍結砂土的屈服強度，對應的橫坐標視為屈服應變，已有試驗結果[25]表明，在相同應變率和負溫條件下，隨著圍壓等級的提高，凍結砂土的屈服強度和峰值應力逐漸增大，而屈服應變、峰值應變和極限應變無明顯變化。因此，在式(6)中引入圍壓增強因子以體現(xiàn)主動圍壓的增強效應：

式中：p為主動圍壓，p0為圍壓對比項，φ為凍結砂土的材料參數(shù)，λ 為圍壓增強因數(shù)。

將式(7)代入式(6)，即可得到考慮圍壓效應的凍結砂土動態(tài)本構模型：

2 凍結砂土動態(tài)本構模型參數(shù)的確定

2.1 確定方法

推導的本構模型共有12 個參數(shù)，其確定順序和方法如下：首先，圍壓對比項p0可根據(jù)試驗條件進行設置，本次設置p0=0.1 MPa；其次，根據(jù)相同負溫和應變率試驗條件下凍結砂土動態(tài)峰值應力與圍壓的關系，可擬合確定圍壓增強系數(shù)λ 和凍結砂土的材料參數(shù)φ；再次，通過將不同應變率條件下的應力-應變曲線相減，能夠得到E2和θ2的范圍；而非線性體的4 個參數(shù)E0、E1、α 和β 主要控制應力-應變曲線彈性段的增長速率，可利用Matlab 軟件中的lsqcurvefit 命令，通過自定義曲線，擬合得到較優(yōu)數(shù)值；然后，塑性體的屈服強度σs主要控制屈服點的位置，其值可根據(jù)試驗得到的應力-應變曲線選取確定；最后，凍結砂土材料的損傷參數(shù)μ和η 需要通過試算，以確定二者對應力-應變曲線的影響，而后通過擬合確定。

2.2 損傷參數(shù) μ 和 η 對應力-應變曲線的影響

采用式(8)的動態(tài)本構模型，以已有的凍結砂土動態(tài)試驗數(shù)據(jù)為基礎[25]，研究了溫度為-15 ℃、應變率為160 s-1、主動圍壓為1.0 MPa 時，損傷參數(shù)μ和η 對凍結砂土動態(tài)應力-應變曲線的影響。由于主動圍壓固定，因此圍壓對比項p0、圍壓增強因數(shù)λ 和材料參數(shù)φ可不設置，模型其余參數(shù)經與試驗數(shù)據(jù)擬合后確定為：E0=5 030 MPa，E1=60 MPa，E2=210 MPa，α=-740 GPa，β=980 GPa，θ2=0.1 ms，σs=9.52 MPa。

根據(jù)上述模型參數(shù)，獲得了不同μ和η 時凍結砂土的動態(tài)應力-應變曲線，如圖3 所示?？梢钥闯觯?1)不同μ和η 條件下，應力-應變曲線的彈性階段和屈服點基本不變，μ和η 主要影響塑性階段和破壞階段曲線的特征。(2)當η=11.5 時，μ值的改變對曲線形狀的影響不大。隨著μ的增大，塑性階段占比逐漸增大，當μ由0.04 依次增大到0.044、0.048 和0.052 時，曲線的峰值應力由11.48 MPa 增大到11.57、11.66 和11.73 MPa，增幅分別為0.78%、1.57%和2.18%，整體增幅較?。欢逯祽儎t由0.026 3 增大到0.029 3、0.031 3 和0.033 3，增幅分別為11.4%、19.0%和26.6%，出現(xiàn)小幅度增大。因此，μ主要影響曲線的峰值應變和塑性段的占比，而對峰值應力的影響較小。(3)當μ=0.048 時，η 的改變對曲線形狀的影響較明顯。隨著η 的增大，曲線中塑性階段占比逐漸增大，且峰后階段即破壞階段的下降速率明顯增大。當η 由5.5 依次增大到9.5、13.5 和17.5 時，峰值應力由11.22 MPa 增大到11.57、11.72 和11.8 MPa，增幅分別為3.12%、4.46%和5.17%；而峰值應變則由0.023 2 增大到0.030 3、0.033 3 和0.035 4，增幅分別為30.6%、43.5%和52.6%，增幅明顯；此外，當η＞13.5 后，η 的變化對曲線的影響逐漸減弱。

圖3 損傷參數(shù)μ和η 對凍結砂土動態(tài)應力-應變曲線的影響Fig.3 Effects of the damage parameters μ and η on the dynamic stress-strain curves of frozen sandy soil

2.3 確定的模型參數(shù)

根據(jù)以上模型參數(shù)的確定方法，并結合損傷參數(shù)對應力-應變曲線的影響，對模型參數(shù)進行試算，最終確定了不同負溫和主動圍壓下凍結砂土的動態(tài)本構模型參數(shù)，見表1。試驗條件為[25]：應變率為160 s-1，溫度為-5 ℃時的主動圍壓分別為0.5、1.0、1.5 和2.0 MPa，溫度為-15 ℃時的主動圍壓分別為0.5、1.0 和2.0 MPa。

表1 凍結砂土動態(tài)本構模型參數(shù)Table 1 Dynamic constitutive model parameters of frozen sandy soil

3 凍結砂土動態(tài)本構模型與試驗結果的對比

3.1 本構模型與試驗應力-應變曲線的對比

圖4 為不同主動圍壓下凍結砂土動態(tài)本構模型曲線和試驗結果[25]的對比，可以看出，不同主動圍壓和負溫下，本構模型曲線和試驗結果[25]吻合較好。本文中推導的本構模型能夠體現(xiàn)圍壓作用下凍結砂土動態(tài)應力-應變曲線的彈性、塑性和破壞3 階段變化特征，也能夠反映曲線中塑性階段占比大以及彈塑性階段拐點明顯的特點，此外，該模型能夠預測主動圍壓對峰值應力的增強作用。

圖4 不同主動圍壓下本構模型曲線與試驗結果的對比Fig.4 Comparison of dynamic stress-strain curves by the modified constitutive model with the test results[25]under different active confining pressures

3.2 本構預測誤差分析

為分析本構模型對試驗結果的預測誤差，以應變?yōu)闄M坐標，以本構模型應力和試驗應力的差為縱坐標，繪制了不同負溫和主動圍壓下凍結砂土的應力差-應變曲線，見圖5?？傮w來說，本構模型應力與試驗應力的差較小，絕大部分可控制在±1 MPa 范圍內，其中最大正誤差出現(xiàn)在彈性階段（對應試驗條件為溫度-15 ℃、主動圍壓1.0 MPa），為1.79 MPa，最大負誤差出現(xiàn)在破壞階段（對應試驗條件為溫度-15 ℃、主動圍壓0.5 MPa），為-1.16 MPa。分析認為，凍結砂土動態(tài)應力-應變曲線彈性段的斜率主要受E0、E1、α 和β 的共同影響，而塑性和破壞階段主要受σs、μ和η 的控制。為進一步分析本構模型和試驗結果的誤差程度，計算了不同試驗條件下凍結砂土的模型應力與試驗應力的平均絕對誤差、均方根誤差和標準差，見表2。溫度為-15 ℃、主動圍壓為2.0 MPa 時的誤差最小，平均絕對誤差、均方根誤差和標準差分別為0.125 MPa、0.193 MPa 和0.158 MPa；平均絕對誤差最大為0.677 MPa（對應試驗條件為溫度-5 ℃、主動圍壓1.5 MPa），均方根誤差最大為1.034 MPa（對應試驗條件為溫度-15 ℃、主動圍壓1.0 MPa），標準差最大為0.794 MPa（對應試驗條件為溫度-15 ℃、主動圍壓1.0 MPa）。

圖5 不同溫度和主動圍壓下凍結砂土應力誤差-應變曲線Fig.5 Stress error-strain curves of frozen sandy soil at different temperatures and active confining pressures

表2 不同試驗條件下凍結砂土的模型應力與試驗應力的平均絕對誤差、均方根誤差和標準差Table 2 Mean absolute errors, root mean square errors and standard deviations between the dynamical stresses of frozen sandy soil by the established model and the test ones under different test conditions

對于動態(tài)應力-應變曲線，峰值應力、屈服強度、峰值應變和屈服應變是凍結砂土最重要的力學和變形參數(shù)，這些參數(shù)對曲線特征的影響顯著。因此，將本構模型對峰值應力、屈服強度、峰值應變和屈服應變的預測結果與試驗數(shù)據(jù)進行對比，見圖6。不同負溫和主動圍壓等級下，本構模型對峰值應力的預測效果最佳，模型峰值應力與試驗峰值應力的最大誤差僅為0.13 MPa；其次預測效果較好的為屈服強度，模型屈服強度與試驗屈服強度的最大誤差為0.6 MPa；模型峰值應變和屈服應變與相應試驗結果的最大誤差分別為0.003 9 和0.002 86，預測效果有待進一步改進。

圖6 不同溫度和主動圍壓下本構模型對動態(tài)參數(shù)的預測結果與試驗結果的對比Fig.6 Comparison of predictions of dynamic parameters by the established constitutive model with the corresponding test results at different temperatures and active confining pressures

4 結 論

基于經典的ZWT 模型，通過引入塑性模型元件，建立了能夠考慮塑性變形的凍結砂土動態(tài)本構模型，研究了損傷參數(shù)μ和η 對動態(tài)曲線和特征點的影響，分析了不同溫度和主動圍壓等級下模型的預測效果，主要結論如下。

(1)本文中建立的本構模型能夠較好地預測主動圍壓作用下凍結砂土動態(tài)應力-應變曲線的強度和變形特性，模型曲線具有明顯的彈性、塑性和峰后破壞3 個階段，此外，模型能夠反映曲線的塑性階段占比大以及彈塑性階段拐點明顯的特點。

(2)損傷參數(shù)μ和η 對動態(tài)應力-應變曲線彈性階段和屈服點的影響很小，主要影響曲線塑性階段和破壞階段的變化，二者的增大都會導致塑性階段占比增大；μ對應力-應變曲線形狀和峰值應力的影響較小，當η=11.5 時，隨μ的增大，峰值應變增大，最大增幅為26.6%；η 對應力-應變曲線形狀和峰值應變的影響較明顯，當μ=0.048 時，隨η 的增大，曲線破壞階段的下降速率明顯加快，峰值應變增幅較明顯，最高可達52.6%，η＞13.5 后，η 的變化對曲線的影響逐漸減弱。

(3)在本次試驗條件下，模型應力與試驗應力的差較小，絕大部分可控制在±1 MPa 范圍內；模型對峰值應力的預測效果最好，試驗峰值應力的最大誤差僅為0.13 MPa；模型對屈服強度、屈服應變和峰值應變的預測效果依次變差。