李宗利,姚希望,李云波,吳正橋,肖帥鵬,劉士達
(1. 西北農林科技大學旱區農業水土工程教育部重點實驗室,楊凌 712100;2. 西北農林科技大學水利與建筑工程學院,楊凌 712100;3. 中水北方勘測設計研究有限責任公司,天津 300222)
彈性地基梁理論因其反映基礎和結構協調變形的優越性,已被諸多學者用于分析寒冷地區土體凍脹和上部結構之間的力學響應[1-5]。Rajani等[1-4]基于彈性地基梁理論分析管道在基礎凍土凍脹時的力學響應。董建華等[5]基于Winkler彈性地基梁理論,研究了寒區錨桿與格構梁復合結構在凍脹作用下力學響應。相比于渠道襯砌凍脹材料力學模型[6-9],彈性地基梁力學模型能反映出襯砌與凍土之間的協調變形[10]。肖旻等[11]基于凍脹力和凍脹強度成線性關系并結合彈性地基梁理論,建立了考慮凍土與結構相互作用的梯形渠道凍脹破壞彈性地基梁模型。李宗利等[12-13]基于自由凍脹量和彈性地基梁理論,結合梯形渠道襯砌凍脹變形特點,給出合理的邊界條件,分別建立了基土均勻自由凍脹和不均勻自由凍脹時的梯形渠道襯砌凍脹彈性地基梁力學模型,計算結果與已有試驗和數值模擬吻合較好。
將彈性地基梁理論引入到凍土凍脹和上部結構相互作用分析中,雖解決了考慮凍土與上部結構之間相互作用的方法問題,但如何合理確定基床系數仍是核心問題。本文用彈性地基梁理論分析凍土凍脹和襯砌板之間的相互作用時,將彈性地基梁理論中的基床系數稱為凍脹反力系數。凍脹反力系數與基床系數存在一定相似性,但卻有所不同:結構對地基施加一定的荷載,地基產生相應的壓縮位移,從而地基給結構一定的反力,地基反力與壓縮位移的比值就是基床系數;基土產生凍脹變形,其上部結構對凍脹變形存在約束,從而形成一定的凍脹反力,凍脹反力和被約束凍脹量的比值即凍脹反力系數。若直接用基床系數分析凍土凍脹和上部結構之間的力學響應有所不妥。同時,考慮到襯砌板下不同點的凍脹反力系數隨被約束凍脹量的不同而發生變化,且每一點凍脹反力系數具體值又是被約束凍脹量的非線性函數,若采用常凍脹反力系數進行計算,計算結果和實際結果會有所偏差。
凍脹反力系數的合理取值影響到基于彈性地基梁理論分析凍脹問題的合理性。目前對基床系數的研究相對較多[14-15],但對凍脹反力系數研究相對較少。同時,基于彈性地基梁理論在渠道襯砌凍脹力學響應上的分析[11-13]尚未考慮到凍土凍脹變形過程中的非線性。本文基于凍土三軸試驗結果,建立考慮圍壓和溫度的鄧肯-張本構模型;參考室內三軸試驗測定基床系數方法,基于數值模擬建立變凍脹反力系數計算式;基于有限差分法離散彈性地基梁平衡微分方程,由此建立變凍脹反力系數梯形渠道襯砌凍脹彈性地基梁力學模型。探究凍脹反力系數分別為變量與常量時在梯形渠道襯砌凍脹力學響應計算結果上的差異,以期為大型梯形渠道襯砌抗凍脹設計提供參考。
1.1.1 基本概念
基床系數是地基土在外力作用下產生單位變形時所需的應力,一般可表示為
式中k為基床系數,MPa/m;P為地基土所受應力,MPa;s為地基變形,m。
1.1.2 三軸試驗確定法
現階段直接通過現場試驗確定地基土的基床系數方法尚不成熟,而室內三軸試驗操作簡單,可控性強。操作過程模擬現場K30原位平板荷載試驗[16]。根據相關規范[16-17],基床系數室內三軸試驗測定法是將土樣經飽和處理后,在有側向靜止土壓力狀態下進行排水固結,側限圍壓應按下式計算:
式中σ1為軸向壓力,MPa;σ3為側向圍壓,MPa;γg為土體容重,kN/m3;hg為土壤所處的深度,m;K0為土體靜止側壓力系數,可按下式計算:
式中ν為泊松比。
固結穩定后,控制圍壓增量Δσ3與主應力增量Δσ1的比值n為某一固定數值,得到Δσ1~Δh0曲線,求得初始切線斜率或某一割線斜率定義為基床系數k。
凍土頂面存在襯砌等結構約束,當凍結達到穩定狀態時,可能出現如圖1所示幾種情況。在上部結構約束作用下自由凍脹量不能被釋放,部分凍脹量被約束,甚至凍土出現壓縮情況,被約束凍脹量的大小直接影響著凍脹反力大小。本文根據被約束凍脹量和自由凍脹量之間的關系將凍土凍脹類型分為自由凍脹、部分約束凍脹、完全約束凍脹和超約束凍脹。
因凍土受力變形具有較強的非線性[18],被約束凍脹量和凍脹反力并不是成正比的關系,如圖2所示。為了建立凍脹反力Pf和被約束凍脹量y之間關系,本文定義凍脹反力Pf和被約束凍脹量y之比為凍脹反力系數kf,即圖2中割線的斜率,如式(4)所示。
式中kf為凍脹反力系數,MPa/m;Pf為凍脹反力,MPa;y為被約束凍脹量,m。
本文定義初始凍脹反力kf0為圖2曲線的初始切線斜率;完全約束凍脹反力系數kf2則定義為圖1中完全約束凍脹情況下凍脹反力Pf2與被約束凍脹量y2(y2=Δh)的比值。
本文參考基床系數室內三軸試驗確定法[16-17],基于數值模擬對凍土的凍脹反力系數進行研究。
凍土凍脹是一個較為復雜的物理現象,涉及到溫度傳遞、水分運移、冰水相變、水結冰后的體積膨脹引起的凍脹。通過數值模擬方法模擬凍土凍脹的全部過程較為困難。因此本文研究的重點是凍土凍結和凍脹達到某種穩定狀態時凍土非線性變形特性對基床系數的影響,暫不考慮凍土凍結和凍脹過程中凍土和襯砌板之間的相互作用。本文以封閉系統飽和凍土為研究對象,忽略水分遷移,認為水分在原位發生完全凍結,溫度對凍土力學性能的影響占主導作用。
1.3.1 溫度傳導方程
凍土凍結和凍脹達到穩定狀態時可按穩態溫度傳導方程對溫度場進行求解。穩態熱傳導的偏微分方程為
式中T為溫度,℃;λ為凍土等效導熱系數,W/(m·℃);?為拉普拉斯算子。
對于凍土凍結和凍脹的穩定狀態,可認為凍土中的水全部凍結為冰。因此對于本文研究的封閉系統飽和凍土可認為是由冰和土顆粒兩相組成。導熱系數等效值可按兩相含量來計算[19],表達式為
式中λi和λs分別為冰和土顆粒的導熱系數,W/(m·℃);θi和θs分別是體積含冰量和體積土顆粒含量,m3/m3。
1.3.2 凍土本構模型
鄧肯-張模型[18]最早是基于非凍土的變形特性提出的,后來一些學者[20-22]將其應用到凍土的變形研究中。本文采用鄧肯-張模型來描述凍土在發生凍脹時的應力應變關系。
式中a和b為試驗常數,εa為軸向應變。
凍土的應力張量S可以分解為球應力張量Tr(εs)和偏應力張量Dev(εs)
式中Tr(εs)為球應變張量;Dev(εs)為偏應變張量;δ為單位張量;K為體積模量,Pa;Gs為割線剪切模量,Pa;εs為應變張量。
鄧肯-張模型的割線剪切模量Gs和體積模量K可按下式計算:
式中G為剪切模量,Pa;γ為剪切應變;qult為極限偏差應力,Pa;E0為初始彈性模量,Pa。
初始彈性模量E0和極限偏差應力qult與鄧肯-張模型中a、b參數有關。當a、b參數與溫度T和圍壓σ3c有關時,初始彈性模量E0和極限偏差應力qult可按下式計算:
式中T為溫度,℃;σ3c為圍壓,Pa。
在考慮凍土凍脹時,應變張量εs可按下式計算:
式中ε為總應變張量;εf為凍脹應變張量。
本文研究的封閉系統的飽和凍土凍脹應變[23]εf為
式中θ為含冰量,m3/m3;ρw為水的密度,ρw=1 000 kg/m3;ρi為冰的密度,ρi=917 kg/m3。
1.3.3 數值模型參數取值及驗證
溫度場較為簡單,僅涉及到冰和土顆粒的導熱系數。冰的導熱系數λi=2.22 W/(m·℃),土顆粒的導熱系數
λs=1.5 W/(m·℃)[19]。
溫度和圍壓對凍土的鄧肯-張模型參數有顯著的影響,對已有試驗數據[22]重新整理,建立考慮溫度和圍壓的凍土鄧肯-張本構模型,再將其帶入COMSOL軟件計算,與試驗結果[22]進行對比,驗證其合理性。
1)凍土三軸試驗
文獻[22]對凍結粉質黏土進行溫度分別為-5、-10和-15 ℃,圍壓分別為0.6、1.0和1.4 MPa的三軸剪切試驗。試樣為Φ50 mm×100 mm的土樣。先將土樣置于真空飽和缸中10 h,使其充分飽和;再對其進行給定圍壓下排水固結;固結完成后進行給定溫度下不少于24 h凍結;凍結結束進行三軸剪切試驗,試驗結果如圖3所示。
2)鄧肯-張模型參數確定及驗證
對該試驗結果進行重新處理,得到不同圍壓和不同溫度下的鄧肯-張模型參數,如表1所示。
對表1中的數據進行回歸分析,可得
式中A1、A2、A3、A4、A5、A6、B1、B2、B3、B4、B5和B6為多項式系數,如表2所示。

表2 鄧肯-張模型a、b參數擬合多項式系數Table 2 a, b parameter fitting polynomial coefficients of Duncan-Chang model
將式(16)~式(17)代入式(9)~式(13),再代入式(8),便可得到考慮圍壓σ3c和溫度T的凍土鄧肯-張本構模型。將其代入COMSOL軟件對該試驗[22]進行數值模擬,結果如圖3所示,可以看出數值模擬結果和試驗結果吻合較好,決定系數R2平均為0.97,說明本文建立的鄧肯-張模型能夠很好地反映文獻[22]試驗凍土樣隨溫度和圍壓變化的力學行為。
1.4.1 三軸試驗數值模型
對于常規三軸試驗,試樣尺寸相對較小,試驗結果只能代表凍土層某一個點的力學特性。本文認為渠道襯砌板所受到的凍脹力全部由凍土層凍脹變形后被約束凍脹量產生。為使測得的凍脹反力系數能綜合反映出整個凍土層的力學行為,數值模擬模型以整個凍土層為研究對象,荷載邊界依然采用規范中室內三軸試驗測基床系數方法中的類似邊界。
如圖4所示,取高度1.5 m、直徑0.5 m的凍土圓柱。土柱為封閉系統飽和土樣,孔隙率約0.45。土柱上表面環境溫度分別為-5、-10和-15 ℃,底部為凍深線所在處,可認為該處為0 ℃。如圖4a所示,數值模擬模型中考慮土柱自重G對應力的影響,在土柱的側面施加K0γgh的圍壓,模擬土柱凍脹前的受力狀態。如圖4b所示,參考基床系數室內三軸試驗測定方法,在土柱凍結凍脹的同時,向土柱頂部施加σ1荷載,在側面施加σ3荷載,逐級增加荷載σ1和σ3,保持增幅Δσ3=mΔσ1,m=0.1[24-25];記錄被約束凍脹量和σ1,此時σ1就是凍脹反力,如圖5所示;根據式(4)和圖2,圖5中曲線的割線斜率即為凍脹反力系數,如圖6所示。根據數值模擬計算結果,土樣在自由凍脹時,垂直方向將會出現61.4 mm的自由凍脹量。
由圖5和圖6可以看出,相同凍結溫度下,被約束凍脹量每增加1 cm,凍脹反力平均增加0.15 MPa,凍脹反力系數平均減小1.12 MPa/m。相同約束凍脹量下,凍結溫度每降低5 ℃,凍脹反力平均增大0.21倍。采用雙曲函數反映這種變化規律,方便工程應用。根據基床系數k取值方法,一般以初始切線斜率或某一割線斜率定義基床系數k[16-17],若以此方法對凍脹反力系數kf取值,明顯會使凍脹反力計算結果偏大。因此,本文認為凍脹反力系數隨被約束凍脹量的改變而發生變化。
1.4.2 凍脹反力系數計算模型
王洪新等[26]認為基坑被動區土體的非線性彈簧地基模型滿足雙曲函數關系,得到不同變形的基床系數計算表達式。該成果對本文整理凍脹反力系數計算結果具有一定的借鑒意義。隨著被約束凍脹量的增加,凍脹反力的增幅逐漸減小,這種規律可通過雙曲線函數反映。因此,本文基于雙曲函數模型對所得凍脹反力和凍脹反力系數結果進行擬合,其擬合公式如式(18)和式(19)所示,擬合參數如表3所示。
式中a′和b′為擬合參數。表3 為本研究土樣所得到的結果,對于其他特定凍土試樣可參考本文的方法進行測定。若條件允許可,通過物理試驗更準確地對凍脹反力系數進行測定。

表3 凍脹反力系數雙曲函數模型擬合參數Table 3 Fitting parameters of hyperbolic function model of frost heave reaction force coefficient
彈性地基梁的基床系數為非常數時,通常很難給出解析解,只能轉而尋找近似的數值解。蔡四維[27]首次將有限差分法用于求解彈性地基梁,用一組差分方程替代微分方程,求解線性代數方程組,使求解過程大大簡化。本節采用有限差分法和牛頓法,通過MATLAB編程,建立變凍脹反力系數的梯形渠道襯砌凍脹彈性地基梁力學模型。
當彈性地基梁上分布荷載q(x) = 0時,彈性地基梁基本微分方程為
式中y為彈性地基梁的撓度,對于本文研究的渠道襯砌凍脹問題即為前述被約束凍脹量,m;E為渠道襯砌板彈性模量,Pa;I為渠道襯砌板截面慣性矩,m4。
如圖7所示,根據有限差分法原理[28],地基梁和其邊界可以被等分成n+4段,共計n+5個節點。設梁撓曲線方程為y=f(x)(被約束凍脹量方程),在撓曲線上取等距離的5個點,間距為Δx,其編號分別為(i-2)、(i-1)、i、(i+1)和(i+2)。
根據有限差分法[28]和材料力學[29]相關原理,節點四階微分和內力的差分方程分別為
式中M為彎矩,N·m;V為剪力,N;Δx為單段長度,m。
將式(18)代入式(20)再代入式(21),可得
式中α=(Δx)4/(EI)。
則圖7中梁上的每一個節點均滿足式(24)。
2.2.1 邊坡襯砌板
根據文獻[12]邊坡襯砌板在凍脹時邊界條件,有
式中Ms為邊坡襯砌板彎矩,N·m;Vs為邊坡襯砌板剪力,N;Fs為邊坡襯砌板在坡腳處受到渠底襯砌板的法向約束力,N;Ls為邊坡襯砌板長度,m。
將式(25)~式(28)代入式(22)~式(23),則有
當x=0(i=0)時
當x=Ls(i=n)時
聯立式(29)和式(30)可得
式中βs=2Fs(Δx)3/(EI)。
聯立式(31)和式(32)可得
2.2.2 渠底襯砌板
根據文獻[12],渠底襯砌板在凍脹時邊界條件為
式中Mb為渠底襯砌板彎矩,N·m;Vb為渠底襯砌板剪力,N;Fb為渠底襯砌板在坡腳處受到邊坡襯砌板的法向約束力,N;Lb為渠底襯砌板長度,m。
和邊坡襯砌板同理,可得
式中βb=2Fb(Δx)3/(EI)。
2.3.1 邊坡襯砌板
將式(33)~(36)帶入式(24)中,圖7中的每一個節點有
式(45)為非線性方程組,可寫成下面矩陣形式
式中M s為邊坡襯砌板系數矩陣;ys為邊坡襯砌板撓度矩陣;βs為邊坡襯砌板邊界矩陣。各矩陣表達式為
2.3.2 渠底襯砌板
和邊坡襯砌板同理,將式(41)~式(44)代入式(24),整理成矩陣形式,則有
式中M b為渠底襯砌板系數矩陣;yb為渠底襯砌板撓度矩陣;βb為渠底襯砌板邊界矩陣。各矩陣表達式為
式(46)和式(50)為非線性方程組,可以通過牛頓法進行計算[30]。牛頓法實際上是一種線性化方法,基本思想是將非線性方程組F(y)=0逐步歸結為線性方程組的求解。
設方程組F(y)=0有近似根ky,將方程組F(y)在ky處展開,有
于是方程組F(y)=0可近似地表示為
這是一個線性方程組,記根為k+1y,則k+1y的計算公式為
將本文式(46)和式(50)代入式(56),則有
式(57)便可采用迭代進行計算,ky和k+1y分別為舊的和新的迭代計算結果;kJ為雅克比矩陣,按下式計算:
對于首次迭代計算
將式(59)代入式(57)進行迭代計算。通過判斷k+1y和ky之差的一次范數是否達到給定精度10-6進行停止迭代,通過MATLAB編程實現這一計算過程。本文數值模型收斂速度快,精度高。如圖8所示,在梁長度1 m、分段數100情況下進行計算,只需經過5次迭代計算,便可完成計算,且具有非常高的計算精度。
文獻[12]根據彈性地基梁短梁解析解,建立梯形渠道襯砌凍脹力學模型。基于該解析解對本文數值解進行驗證。由于彎矩和凍脹反力可根據撓度(凍脹量)進行求解,因此只需對撓度(凍脹量)進行驗證即可。對文獻[12]中的工程算例采用解析解和數值解進行分別計算,結果如圖9所示。可以看出本文基于數值模型計算出的結果和解析解基本吻合。
襯砌板長度與厚度比值體現了襯砌的抗彎剛度,反映了襯砌板對渠基凍土凍脹約束能力大小,影響凍脹量分布,而渠道襯砌板厚度取值范圍不大,因此結合襯砌長度來分析凍土非線性變形特性影響規律,即襯砌剛度對渠道襯砌凍脹的影響。
本節以某寒區梯形襯砌渠道為例,探究凍脹反力系數對渠道襯砌凍脹的影響。渠道基土自由凍脹量5 cm,環境溫度-15 ℃,襯砌板厚度10 cm,襯砌混凝土彈性模量24 GPa。凍脹反力系數為變量時,選取本文-15 ℃計算結果;凍脹反力系數為常量時,選取-15 ℃計算結果初始值。邊坡襯砌板和渠底襯砌板計算結果如圖10和圖11所示。
3.1.1 凍脹量
如圖10a所示,常凍脹反力系數凍脹量計算結果比變凍脹反力系數凍脹量計算結果略大,但分布規律基本一致。當邊坡襯砌板長度為1 m時,凍脹量分布接近直線,隨著長度增加,凍脹量曲率逐漸增加。這主要是因為在自由凍脹量一定,邊坡襯砌板坡腳約束不變的情況下,隨著邊坡襯砌板長度的增加其整體剛度逐漸變小。
3.1.2 凍脹反力
如圖10b所示,凍脹反力系數分別為常量和變量時,計算結果在法向凍結力的分布和大小上基本相同,但在法向凍脹力上卻有較大差異。變凍脹反力系數計算出的法向凍脹力明顯比常凍脹反力結果小,同時距離坡腳越近,被約束凍脹量越大,計算結果差異越明顯。本文對梯形渠道襯砌凍脹過程簡化,認為在渠底襯砌板橫向約束作用下,邊坡襯砌板坡腳處的自由凍脹量完全被約束。自由凍脹量相同時,不同長度邊坡襯砌板在坡腳處被約束凍脹量相同。由Winkler地基理論[28]可知,若地基梁上某點變形相同,則產生的反力相同。但由于凍脹反力系數的差異,常量凍脹反力系數計算出的最大凍脹反力是變量的1.43倍。
3.1.3 彎矩
由圖10c可以看出變凍脹反力系數彎矩計算結果明顯小于常凍脹反力系數計算結果;變凍脹反力系數計算彎矩最大值位置相對于常凍脹反力系數整體略遠離坡腳。隨著長度增加,彎矩最大值逐漸變大后趨于穩定。這主要是因為在渠基凍土在自由凍脹量一定時,隨著邊坡襯砌板長度增加,整體剛度逐漸變小,遠離坡腳區域的襯砌板對渠基凍土凍脹約束較小,隨基土凍脹一起發生上抬位移;在坡腳處由于渠底襯砌板對邊坡襯砌板坡腳處的約束,使得該處產生較大的凍脹反力,但邊坡襯砌板只有坡腳點處的凍脹受到來自渠底襯砌板的強迫約束,其影響范圍相對有限,且當邊坡襯砌板長度增加到一定程度時,影響范圍基本不變;在被約束凍脹量一定的情況下,影響范圍內產生的法向凍脹力也基本相同,因此當邊坡襯砌板長度增加到一定值時,其彎矩最大值基本保持不變。常量凍脹反力系數計算出的彎矩最大值平均是變量的1.12倍。
3.2.1 凍脹量
如圖11a所示,在自由凍脹量一定的情況下,隨著襯砌板長度的增大,中間部位被約束凍脹量逐漸減小,變凍脹反力系數和常凍脹反力系數的差異逐漸減小,導致隨著渠底襯砌板長度的增加,變凍脹反力系數和常凍脹反力系數的凍脹量計算結果差異逐漸變小;中間部位凍脹量逐漸增大,但增幅逐漸變小。
3.2.2 凍脹反力
如圖11b所示,當渠底襯砌板長度較小時,變凍脹反力系數的凍脹反力計算結果和常凍脹反力系數的凍脹反力計算結果相差較大,反之相差較小。由于在坡腳被約束凍脹量都是5 cm,因此不同長度的渠底襯砌板計算的凍脹反力在坡腳處相同,和邊坡襯砌板原因相同。結合圖11a在渠底襯砌板兩邊自由凍脹量被約束,中間渠道基土凍脹發生了一定的變形,釋放了一定的凍脹力,因此渠底襯砌板凍脹反力分布為兩邊大,中間小。當渠底襯砌板長度為3和4 m時,在渠底襯砌板中部局部產生了不大于0.15 MPa的法向凍結力,這主要是由于襯砌板長度相對較長時,中間區域產生的凍脹量相對較大,使得襯砌板中間有遠離基土的趨勢,但渠道基土和襯砌板凍結在一起,因此,在中間區域產生了一定的法向凍結力。
3.2.3 彎矩
由圖11c可以看出隨著渠底襯砌板長度的增加,變凍脹反力系數的彎矩計算結果和常凍脹反力系數的彎矩計算結果差異逐漸變小,和凍脹反力出現這種現象的原因基本一致。當渠底襯砌板長度增加到一定程度時,在渠底襯砌板中部出現法向凍結力,如圖11b所示,對中部彎矩有一定削弱作用。因此,當渠底襯砌板長度超過一定值時,隨著長度增加,中部彎矩開始減小,彎矩最大值出現在兩邊。
1)渠道基土在凍脹過程中具有非線性變形特性,基于室內三軸試驗確定基床系數方法是可行的。提出了雙曲函數反映基床系數隨約束凍脹位移的變化規律,方便工程應用。
2)探究了凍脹反力系數分別為變量與常量時,在梯形渠道襯砌凍脹力學響應計算結果上的差異。結果表明,對于邊坡和渠底襯砌板,常量凍脹反力系數計算出的最大凍脹反力是變量的1.43倍,計算出的彎矩最大值平均是變量的1.12倍。因此在采用彈性地基梁理論分析渠道襯砌凍脹問題時,若凍脹反力系數采用為常量,不考慮凍土的非線性變形,會使得計算結果偏大。
3)本文所建立的凍脹反力系數計算數值模型未能很好考慮水分遷移,得到凍脹反力系數比實際值偏小,但不影響所建立的梯形渠道襯砌凍脹彈性地基梁力學模型的應用。同時,尚未考慮襯砌板與地基分離和襯砌板的凍縮對彎曲影響,計算結果和實際值存在偏差,未來模型仍需進一步完善。