熱-濕-應力耦合作用下框架錨桿支撐的邊坡凍融特性分析

劉銘, 鮑碩超

(吉林建筑大學土木工程學院, 長春 130000)

凍土區的存在及其演變過程誘發的熱融滑坡[1]對自然環境和人類生產活動造成了不同程度的影響,直接威脅人民的生命財產安全,與邊坡穩定性密切相關的是邊坡防護工程。傳統的支護結構在凍土邊坡治理中存在諸多問題,如部分防護結構在凍脹作用下破壞,無法達到理想的防護效果。因此,對凍土邊坡防護需要一種安全有效、技術先進、施工方便、造價經濟的方法[2]。

對于凍土區邊坡防治問題,許多學者在凍土邊坡的理論與防治措施方面進行了大量研究。Zhang等[3]研究發現,凍融破壞是凍土區邊坡失穩的最主要原因,其本質是土壤水分在土壤孔隙中的凍結和融化。Xin等[4]針對不同土質進行試驗發現,不同性質的土體經歷凍融循環后其強度變化也不同,土體黏聚力和摩擦角隨著干密度、凍融循環和飽和度的增加而逐漸減少。Abdolghanizadeh等[5]提出凍融循環次數增加導致土體產生的非線性應力松弛,是邊坡穩定性研究的重要部分。Vahdani等[6]在對凍融數值模型的研究中,將溫度作為一個額外的變量,引入一個廣義的應變和應力張量,以此來計算對邊坡穩定性的影響。

近年來,學者們針對邊坡的支護展開了大量研究。何江飛等[7]針對黃土高邊坡治理加固中的空間受限和高填方填筑等問題,提出了加筋土-框錨組合體系,對該組合體系協同作用進行了分析。曹程明等[8]以特殊土巖組合地層為背景,采用樁錨支護,結合監測數據對支護結構內力變化規律進行研究。在眾多支護結構中,框錨結構因其簡單可靠的設計原理和低廉的造價得到了廣泛的發展與應用[9-10]。基于此,以框架錨桿結構支護邊坡,設計室外試驗。在理論上建立凍融荷載作用下邊坡土體的物理場數值計算模型,土體和結構之間工作作用的計算模型。基于模型方程,使用MATLAB編制計算程序,通過與試驗結果的對比,驗證程序的正確性。得到季凍區框錨支護邊坡的水分場、溫度場、應力場以及錨桿內力分布規律,揭示支護結構和土體相互作用的機理,為框架錨桿結構支護凍土區邊坡的應用提供可靠的理論依據。

1 數學模型

1.1 基本假定

(1)土體土質均勻聯系,為各向同性彈性體。

(2)水汽蒸發的耗熱量忽略不計,且不計其他化學勢和場的作用,如鹽分的析出。

(3)不考慮施工時對支護結構造成的誤差,假定其一直處于正常工作情況。

1.2 基本方程

根據多孔材料傳熱傳質理論[11],土體凍融過程中溫度與水分的變化歸結為以下二維方程組[式(1)～式(3)]。

(1)

(2)

ψ=P+G

(3)

式中:C為土體密度ρ和土體比熱容cp乘積,即ρcp;t為計算時間;λ為土體的導熱系數;θi為體積含冰率;θu為體積未凍水率;ρi和ρw分別為冰和水的密度;kw為與未凍水相關的滲透系數;ψ為土水勢能,水壓力P和重力勢G之和。

水壓力P可通過Clapeyron方程[12]描述為

(4)

式(4)中:T0為土體凍結溫度;L為冰水相變潛熱。

由平衡微分方程和幾何方程可得

(5)

式(5)中:?為拉普拉斯算子;σ為應力矩陣;G為體積力;ε為總應變矩陣;u為位移矩陣。

在凍結過程中,土體體積變化是由原土體中部分水和遷移來的部分水凍結成冰引起的。在凍結和融化過程中假定土體顆粒不可壓縮,體膨脹變形εv可表示為

εv=0.09(θ0-θu)+1.09Δθ-n

(6)

σ=c(ε-ε0)

(7)

1.3 結構與土體協調作用

將框架柱等效為梁單元,錨桿等效為桿單元,土體單元采用三角形單元。為保證結構精細化模擬的有效性,至關重要的是單元形成的數值模型的準確與否。因此,采用零厚度四節點單元作為結構和土體的聯結單元[13]。如圖1所示,聯結單元寬度w=0,節點1和4、3和2在開始時占有同一空間位置,可以很方便地放置于結構和土體之間,不影響結構和土體的單元劃分。

1、2、3、4為聯結單元的4個節點序號

聯結單元中上層節點1、2與土體單元聯結,將土體單元的位移和應力狀態傳遞給下層節點3、4,下層節點再傳遞給結構單元。通過聯結單元的形狀函數來描述結構和土體的運動關系,進而體現出二者的聯結關系。

1.4 結構與土體協調作用方程

土體為三角形單元,其位移模式[14]為

(8)

式(8)中:ui、uj、um、vi、vj和vm分別為三角單元3個節點i、j、m上x、y方向的位移;NS為土體單元的形函數矩陣;δeS為土體單元的節點位移。

由虛功原理可得到土體單元的剛度矩陣和節點力與節點位移的關系可表示為

(9)

feS=kSδeS

(10)

(11)

(12)

(13)

對聯結單元其位移可表示為

=NBδeB

(14)

式(14)中:ΔU、ΔV為聯結單元在x、y方向上的位移;Uup、Udown、Vup和Vdown分布為單元上下邊界的在x、y方向上的位移;NB為聯結單元的形函數矩陣;δeB為聯結單元的節點位移。

由虛功原理可得到聯結單元的剛度矩陣和節點力與節點位移的關系可表示為

(15)

feB=kBδeB

(16)

式中:Ω為單元計算空間;kB為聯結單元的剛度矩陣;k′為聯結單元的剛度系數;feB為聯結單元的節點力。

通過對應力場求解得出δeS和feS,即求得聯結單元上層節點位移和節點力,代入式(16)即可求得下層節點位移和節點力即δeB,即與結構單元相對應的節點位移,將其代入結構單元中進行計算即可求得結構的節點力。

feStr=kStrδeB

(17)

式(17)中:feStr為結構單元的節點力;kStr為結構單元的剛度矩陣。

將各部分的單元剛度矩陣在整體坐標系下進行迭加可得到結構與土體的整體協調作用方程為

(18)

(19)

(20)

FS=KSδS

(21)

FB=KBδB

(22)

δB=δStr

(23)

FStr=KStrδStr

(24)

由式(1)、式(2)和式(18)～式(24)組成了季凍區框架錨桿支護邊坡的水熱力耦合模型。

2 建立數值模型

對水分場和溫度場采用有限差分法計算得到每個時間步長的結果,對應力場在每個時間步長內采用有限元法進行分析。前文已經給出應力場的有限元方程,不再贅述。僅對溫度場的差分公式進行離散分析,在對水分場離散時亦是同理。

2.1 有限差分法

有限差分法[15]是以差分的原理為基礎的計算方法,通過離散的函數值所構成的差商來近似逼近相應的偏導數。如溫度T對時間t的導數可表示為

(25)

式(25)中:Tn+1為n+1時刻的溫度;Tn為n時刻的溫度;tn+1和tn分別為n+1時刻和n時刻;Δt為計算的時間步長。

在幾何意義上表示求弧線的斜率,如圖2所示。

圖2 斜率的表示Fig.2 Representation of slope

2.2 建立差分方程

將式(2)～式(4)代入式(1)整理得

(26)

(27)

(28)

為方便理解,在此引入半線。即在點(i,j)和點(i+1,j+1)之間引入點(i+1/2,j)和(i,j+1/2)。式(26)左端中在點(i,j)處T對t的一階偏導等于點(i,j)在n+1時刻時T和n時刻T的差除以時間步長Δt。即

(29)

(30)

式(26)右端中在n時刻點(i+1/2,j)處T對x的一階偏導等于n時刻時點(i+1,j)T值和點(i,j)T的差除以空間步長Δx。即

(31)

(32)

n時刻點(i,j)處T對x的二階偏導等于n時刻時點(i+1/2,j)一階偏導和點(i-1/2,j)一階偏導的差除步長Δx。即

(33)

T對y的二階偏導同理展開即可。

假設在n時刻時,點(i,j)的溫度T已知。對式(26)展開后左端有兩個時間層n+1、n,右端若都以n時刻展開計算求得的n+1時刻結果誤差極大,而都以n+1時刻計算對計算機內存需求增大,變得煩瑣。因此,采用交替計算方法即能解決誤差又能使計算過程簡化。

在n時刻,各點溫度值已知,如圖3所示。

xi-1、xi、xi+1、yi-1、yi、yi+1為n時刻或n+1時刻時,網格單元劃分后的坐標

(34)

對于邊界條件,在上邊界是T關于時間t的函數,即T=T(t);在下邊界T為一個已知值。左右邊界的熱流密度為0,即

(35)

由左邊界條件可知:

(36)

對左邊界點展開,即i=1時。如圖4所示,得到左邊界點的差分方程,右邊界點同理展開即可。

圖4 左邊界節點Fig.4 Left boundary node

(37)

對所有離散點都進行如上展開,即在n+1時刻在x方向上展開,在n時刻在y方向上展開。而n時刻時各點T已知,結合邊界條件寫成矩陣形式為

AU=B

(38)

(39)

(40)

(41)

(42)

(43)

對矩陣方程求解可得到n+1時刻沿x方向展開的離散點T值。

與計算n+1時刻類似,在n+2時刻,此時n+1時刻各點T值已求得。對內部點在y方向上以n+2時刻展開,在x方向上以n+1時刻展開,即能得到n+2時刻的差分方程。不同的是在n+2時刻對邊界的處理。

由上邊界條件可知,在n+2時刻時上邊界點的溫度已知,即j=1時。如圖5所示,得到上邊界點的差分方程,下邊界點同理。

圖5 上邊界節點Fig.5 Upper boundary node

(44)

在n+2時刻,對離散點都進行相應的展開后組合成矩陣,求得n+2時刻的解。如此反復,便可得到所有時間步長內的解。

3 程序驗證

為了得到真實環境下各物理場變化關系,結合室外環境設計了室外試驗,將其結果與計算結果進行對比驗證。

3.1 室外試驗

試驗模型坡高3.0 m,坡角為75°,如圖6所示。框架混凝土強度等級為C30,橫梁尺寸為20 mm(長)×20 mm(寬)×3.0 m(高),立柱尺寸為20 mm(長)×20 mm(寬)×2.0 m(高)。錨桿鋼筋為HRB400級,錨固段注漿采用M20級水泥漿,半徑40 mm。試驗相關設計參數和傳感器布置如圖7、圖8所示。

圖6 框架錨桿支護邊坡室外模型試驗Fig.6 Outdoor model test of frame anchor support slope

圖7 模型示意圖Fig.7 Model sketch map

3.2 參數取值

由地勘報告可知,試驗場地土體為常見的高原地區黃土,土質均勻、松散、欠固結。通過室內試驗得到地基土平均含水率為20.58%,內摩擦角為15°,黏聚力24.87 kPa,孔隙比為0.98。經判斷,地基土屬于弱凍土層,凍結溫度Tf為-0.21 ℃。結合相關資料[16-18],基本特性參數如表1所示。

表1 基本參數Table 1 Basic parameters

凍土導熱系數λf和融土導熱系數λu分別表示為

(45)

凍土熱容Cf和融土熱容Cu分別為

(46)

導水系數kw為

(47)

未凍水與溫度的關系可表示為

(48)

土體彈性模量E與泊松v分別表示為

(49)

(50)

3.3 計算模型

3.3.1 幾何尺寸

結合試驗模型和場地因素考慮,幾何計算模型如圖9所示。其坡高3.0 m,坡角為75°,長為6.8 m,寬為9 m。

圖9 模型尺寸Fig.9 Model size

3.3.2 初始條件

對于溫度邊界條件,由于試驗場地的日平均氣溫會低至為-4.1 ℃,結合室外試驗實測結果,坡頂、坡面和坡腳的熱邊界條件為

(51)

式(51)中:t為時間,d。

土體初始溫度由觀測數據取夏季凍融周期外的平均溫度4.5 ℃。

對于水分邊界條件,由于正午溫度的升高,地表水分蒸導致上邊界水分散失,而降雨降雪導致上邊界水分補充,左右邊界及下邊界無水分滲入或滲出,因此上邊界的水分條件為

(52)

土體初始含水率由土樣室內試驗測定的20.58%。

位移邊界條件:左右邊界在x方向上固定y方向上自由,下邊界x、y方向均為固定,其余邊界自由。

3.4 試驗結果驗證

凍結深度隨時間變化的曲線如圖10所示。可以看出,前半段計算值與試驗值較吻合,后半段與試驗值相比偏大。從凍結深度曲線變化情況看,計算結果符合凍脹隨時間發展的規律,因此本文建立的模型是宜用于季凍區邊坡工程的設計和計算。

圖10 不同時刻下凍結深度Fig.10 Frezzing depths at different time

4 結果與分析

由試驗觀測結果可知,最大凍結深度與最大融化深度分別出現在1月5日和4月4日。其中,凍結時表示凍深達到最大值,融化時表示融深達到最大值。

4.1 溫度場

凍結時刻前后12 d的觀測數據與計算值對比曲線如圖11所示。凍結和融化時的溫度云圖如圖12所示。

圖12 溫度云圖Fig.12 Temperature nephogram

由圖11可知,計算值整體比試驗值小約2 ℃,但試驗值和計算值呈現出相同的發展規律,隨著時間的進行,溫度逐漸降低,凍結深度不斷往下發展。在最大凍深后,溫度有小幅度的回升,可見,其變化符合凍脹規律,計算值與試驗值基本吻合。

由溫度云圖(圖12)可知,凍結最大影響深度約2.0 m,再往下溫度基本保持穩定。受坡頂凍結的影響,坡面頂部凍深最大。融化時地表溫度影響深度約1.5 m,受坡頂融化的影響,坡面融化深度沿坡面減小。結合凍結和融化的溫度云圖可知,在整個凍融過程中,坡體內部只有一定深度的土體溫度發生了變化,這一部分土體經歷凍脹和融化變形,其應力狀態也隨之變化,容易發生凍融滑塌,在實際工程中應重點觀測該部分土體。

4.2 水分場

凍結時刻前后12 d的觀測數據與計算值對比曲線如圖13所示。凍結和融化時的水分云圖如圖14所示。

圖13 未凍水對比Fig.13 Comparison of unfrozen moisture

圖14 未凍水云圖Fig.14 Unfrozen moisture nephogarm

從圖13可以看出,計算值比試驗值要大約0.2%,但與試驗值吻合相當好。隨著凍脹變形的發展,未凍水含量不斷減少直至水分完全凍結,這也符合實際情況。

未凍水含量與溫度相關,由圖14可知,凍結時未凍水云圖與溫度云圖趨勢接近。在2.0 m處分界明顯。在融化土層內,受坡頂融化影響,坡頂未凍水量最低,沿坡面向下水分不斷聚集。

4.3 應力場

凍結與融化時的剪應力云圖如圖15所示。凍結時土體性彈性模量的增加,導致了剪應力是融化時的2.5倍,凍結時剪應力沿坡面向下大致呈線性分布,在坡腳集中。在邊坡其他區域剪應力分布均勻,因此在凍結時邊坡整體是較穩定狀態;融化時,剪應力在坡腳出現了明顯的應力集中,在融化區域和未融化區域界面發生剪應力突變,邊坡處于不穩定狀態。

4.4 錨桿軸力分析

錨桿軸力如圖16所示。在初始、凍結、融化3種工況下,結合試驗結果與計算值對比分析可知,相比于初始時,由于在邊坡坡面產生的凍脹變形,進而導致支護結構發生位移,直觀的表現在第一排錨桿自由端軸力增加了約20 kN,第二排錨桿自由端軸力增加了約17 kN,并沿著錨桿長度方向變化量不斷減小。離坡面越遠,這種由凍脹變形引起結構的位移量越小,因此在錨固段軸力變化不大。3種工況下,凍結時錨桿軸力最大,融化后軸力大幅度降低,但仍然要比初始時的軸力值要大,而造成這種現象的原因是融化后土體應力發生了重分布,坡面產生了不可恢復的位移,從而導致錨桿無法恢復至初始工況下的軸力。

圖16 錨桿軸力對比Fig.16 Comparison of Anchor Rod Axial Force

5 結論

(1)建立框架錨桿支護邊坡的耦合計算模型,結合MATLAB編制了計算程序,并與室外試驗對比分析,驗證了程序的可靠性。

(2)在試驗過程中存在只有一定深度的土體僅受到凍融影響,這一深度可以通過計算得到。與凍結時相比,融化時錨桿軸力大幅度降低,但仍留有殘余應力,表明凍融作用對結構內力影響較大,在采用結構支護邊坡時,應按凍脹工況對結構進行設計。

(3)考慮多物理場耦合作用,通過試驗和理論計算,對框錨支護邊坡在單次凍融條件下各物理場反映的規律進行分析,但并沒有進一步分析多次凍融循環條件下各物理場的變化,需在以后的工作中對此部分進行深入研究。