1、凍土損傷本構模型和耦合問題的數值模擬作者：朱志武， 寧建國，馬巍作者單位： 西南交通大學力學與工程學院，成都610031 ;北京理工大學爆炸科學與技術國家重點實驗室， 北京 100081 ；中國科學院寒區旱區環境與工程研究所凍土工程國家重點實驗室，蘭州730000收稿日期： 2010-01-06； 接受日期： 2010-02-23摘要從材料的微觀力學機理出發， 建立了凍土損傷的彈性本構模型。 對于不同溫度和冰體積 含量下的凍結砂土， 由該模型計算的結果與實測的應力應變曲線較吻合， 建立的凍土損傷 本構模型能夠較好地描述實際凍土材料的力學性質。 利用自行開發的有限元程序，加入所 推出的本構關系和

2、所建立的數學力學模型，通過對渠道凍結和凍土路基的水分場、溫度場、 應力場進行數值模擬計算，得到了比較準確、詳盡的符合實際的溫度場與應力場、位移場、 應變場耦合的計算結果，與前人的計算及實測結果相吻合，規律一致。算例表明該程序能夠計算凍土材料的相關物理量，且能很好地描述它們之間的關系。研究在前人的工作基礎上，結合青藏鐵路路基工程， 通過對凍土水、熱、力三場耦合機理這一固體力學學科領域 的難點問題研究， 為凍土工程的設計、施工和維護提供有效的科學依據、分析模型及必要 的參考。關鍵詞 凍土， 微觀力學， 路基， 耦合， 應力場1 引言凍土是一種對溫度敏感和易變的特殊低溫地質體， 多年凍土的面積約占全

3、球陸地面積 的 23%， 其中我國多年凍土地區約占國土面積的21 。 5%。 隨著社會的發展，凍土地區的道路工程、水利工程、建筑工程、礦山工程、能源工程等應運而生并迅猛發展，特別是隨著 青藏鐵路的建成， 關于凍土的研究顯得越來越重要而急迫。凍脹融沉最本質的要素是土體 凍融過程中水、 熱、 力三場耦合問題， 但由于問題的復雜性和研究手段的限制， 目前研究 僅以水、熱兩場耦合問題研究的最多， 忽略應力場的問題，這樣會導致凍融過程研究的不 完整， 進而影響到凍脹融沉預測的準確性。 因此，水，熱、力三場耦合問題是凍土研究領 域的關鍵問題所在， 也是國際前沿研究課題， 研究結果可用于凍土構筑物的設計和穩

4、定、 安 全評估中，促進凍土材料的大力發展，具有廣闊的工程應用前景。目前人們對凍土中冰-土的相互作用已經有了比較深刻的認識和研究。 Fremond1 對正凍土的熱力學性質進行了初 步描述；Gary2對凍融循環進行了系統的分析；Konard3提出了一個描述正凍土的冰晶形成與水分遷移的模型；Denn is4對凍結過程進行了深入系統的試驗研究；我國學安維東等人 5先后對凍土的水分遷移與熱質遷移， 水熱力耦合及其本構問題，進行了較深入的試 與理論研究。 然而以往的研究大多從熱力學、混合物理論等角度出發建立起凍土的各種力 學模型， 多重于兩場 （溫度場與水分場 ）的耦合作用， 力場只是在分析凍結溫度時，

5、作為一項計算相變溫度的指標被引進，專門從力學機理出發考慮凍土的三場耦合及其本構關系的研究尚未見報道。而在寒區的工程建設，人們不可避免會遇到諸如凍脹融沉等凍土問題并開展相應的研究工作。Foriero等人用有限元模擬了凍土邊坡的蠕變；Su等人7對多年凍土區通風路基進行了數值分析研究；Dempsey8建立了非飽和土中水熱耦合運移的數學預報模型；Rai等人9導出了流體、熱流和土體變形完全耦合的控制微分方程；Masters10把位移、孔隙壓力和濕度值當作初始未知量，研究了土體變形、流體溫度和孔隙壓力之間的耦合規律；Zhang等人11對寒區碎石路基在凍結過程中邊界條件的影響進行了深入研 究。對寒區工程中地

6、基基礎進行溫度、水分、應力等的分析具有十分重要的現實意義，可以為實際施工和道路的正常運營提供有益的參考。因此，本文針對凍土工程中急待解決的土體凍融過程中水、熱、力三場耦合的力學機理問題開展研究，結合青藏鐵路路基工程，從材料細觀力學出發，建立含損傷的凍土本構模型，在對凍土本構模型的研究中走出了一條新路。根據傳熱學，滲流理論和凍土力，建立了凍土溫度場、水分場、應力場耦合問題的數學力學模型，并對三場之間的耦合作用進行相應的數值模擬研究。2凍土的損傷本構模型2。1凍土的彈性模量從細觀力學的角度出發，首先將凍土看作由土和冰組成的復合體單元。其次，把在整個凍土體中所占比例很大的土顆粒作為骨架，而把冰作為填

7、充體看待，根據復合材料理論中的經典混合律思想將其耦合為凍土的本構關系，最后考慮加入損傷的影響。假設土顆粒為均勻連續體，土顆粒與冰之間完全黏結。令土和冰的彈性模量分別為Es和Ei,剪切模量為Gs和Gi,泊松比為v s和v io下面各式中的下標 s和i分別表示凍土中土和冰的各 個分量。于是由細觀力學的混合律理論可得到凍土復合材料的等效彈性模量K和等效剪切模量G為K(1)式中：和為土和冰的體積模量；和為土和冰的剪切模量。根據各向同性材 料彈性常數之間的關系：K =E 3(1- 2 v ), G = E 2(1+v )。經推導后可得到由土、冰的彈性模量和泊松比表示的凍土等效彈性模量 的表達式為12E

8、=乙& （1-2叫）+q爲（1 一 2匕腑遲（1+叫）+柘（1 +1 瓦（1 +片）（1 2片）兀匚（1芮E和等效泊松比vV(3) 咔憶(1+叫)(1-2叫)剛(1+叫)(1-2叫) G 瓦(1+片)(12h)+G 爲(1 +叫)(1一2叫廠 1C , C式中： 和分別為土和冰的體積含量（也稱為體積分數）。2。2凍土的損傷本構關系根據Lemaitre等效應力原理13，即應力（T作用在受損材料上引起的應變與等效應力d作用在無損材料上引起的應變等價。受損凍土在一維狀態下的應力-應變關系為crE(l D)(5)上式中：d為損傷力學中的等效應力，此等效應力相當于土力學中的有效應力;E為無損材料的彈性模

9、量，即初始彈性模量，由（3）式確定；E為受損材料的彈性模量,即有效彈性模量；D為損傷量。因此，根據損傷力學理論14,得到凍土材料內部損傷型本構關系為 d = (1- D)Es 2。3損傷演化規律在外載作用過程中，由于材料強度服從概率統計上的Weibull分布，因而可以認為材料的損傷量D朱志武等：基于損傷的凍土本構型及水、熱、力三場耦合數值模擬研究也服從該分布。在此借助混凝土損傷的研究結果15，用雙參數的Weibull分布表示凍土的損傷量：經推導，可得式中：e為應變，n和a分別為形狀參數和尺度參數，均為非負數。(8)上式即為以應變作為損傷演化控制量的損傷演化方程， f 為應力峰值所對應的應變值，

10、 n為表征材料損傷演化特征的材料參數。 將（8）式代入（6）式中，從而得到凍土的損傷本構 關系為可以看出，凍土中任意點的應力 d與凍土的彈性模量、極限強度、應變峰值及該點的應變 有關。 凍土的損傷形狀參數 n是彈性模量E和割線模量Em的函數， 一般可由實驗或細 觀力學方法加以確定。從理論上講，n值的大小與凍土材料的性態及外表破壞模式有一定 的關系，即：n值越大，凍土越趨向于脆性破壞；n值越小，凍土越趨向于塑性破壞。2。4計算實例及分析根據已知土和冰的彈性模量，選取冰不同的體積含量，利用（3）和式得到了凍結砂土在不同溫度下的等效彈性模量如表1,凍結砂土的孔隙比 e為0。40。5、天然含水量（w%

11、）為1518。利用表1中的彈性模量值，由（9）式計算得到的應力-應變曲線， 與在相同溫度下單軸壓縮試驗得到的應力-應變曲線16的比較， 如圖1和2所示。表1不同溫度下揀結砂土的零效模険冰體積含貳（）砂七彈性沌沽磁匕彈性一4910U0546005-S 46342.046.0363.2B5-8.06820.0387.221-2.U4910.0534.9570.1-5.404X04.0678 1376820.0726.053一工。4910.01021.25&0.2-5 46342.0斗&01307.694-8.06820.0I40J 3C3凍結砂土在不同溫度下的應力-應變曲線比較(k;=0.05)2

12、考慮損傷后在不同冰體積含量E和不同溫度下計算得到的應力應變曲線G備=005; JK/=0.1 * LK/=0.2(12)可以看出，溫度是影響凍土變形的重要因素。隨著土溫的降低， 土中未凍水含量減少，凍土中冰含量增加，固體顆粒間的膠結力增加，整個凍土體的變形式也由流變和塑性破壞逐漸轉向脆性破壞。在相同的應力水平下，凍結砂土的變形隨著土溫的降低而減小，彈性變形在總變形中所占的比例卻隨著土溫的降低而增大，所占的殘余變形減小。另外， 從圖中還可以知道， 彈性變形盡管隨著凍土溫度的降低而有很大增長，但其所占的比例仍然很小， 特別是在凍土溫度較高的情況下（一般高于-5 C范圍內）,彈性變形只占總變形量10

13、%25%。在凍土溫度較低時，彈性變形所占比例也不會有很大提高。為了更清楚地體現冰的體積百分比i c對凍土損傷強度的影響，我們將（9）式得到的在不同冰體積含量的應力-應變曲線進行比較，如圖2所示。 可以看出，在不同溫度下，由于凍土中冰含量的增加，凍土體的損傷強度也相應提高。但在相同溫度下，隨著應變的增大，凍土體的壓縮應力也逐漸達到最高點（屈服點）并已經開始下降，并且不同冰體積含量的凍土的應力屈服點也非常接近，說明在凍土材料破壞（屈服）時冰的體積含量對凍土體的屈服強度影響不大。3 土體凍結水、熱、力三場耦合數值模擬土體在凍結或融化過程中，溫度、水分、應力三場的相互作用是一個極其復雜的熱力學、物理化

14、學和力學的綜合問題。由于水-冰之間的相變以及水分遷移聚冰現象的存在，土體在凍融過程中溫度場和水份場的變化將伴隨著應力場的變化，同時溫度場的變化也會引起水分遷移和應力狀態的改變；反過來，應力場的變化又引起溫度場和水分場的重分布。這三個場相互制約的關系存在于整個土體的凍融過程中，成為凍土材料本身所特有的水、熱、力三場耦合問題。以自行開發的有限元程序為平臺，加入本文推導出的本構模型，對實際工程中凍土的水分場、溫度場及應力場之間的相互耦合問題進行數值模擬，與工程實測和前人所做的實驗結果進行對比分析。3。1渠道凍結三場耦合數值分析 3。1。1控制微分方程,dT0 +一dt dx1務丿dx*1即丿(10)

15、其值分別為:(11)正凍結過程中渠道基土的水熱耦合，可表述為二維非線性變系數拋物型偏微分方程:式中，c和入為路基土的等效容積熱容量和等效導熱系數,C1-dWucT(13)未凍水含量u W與溫度T的關系由室內試驗確定5:氏=f(nC+,TTbTTp.(14)C-,TTb.A+,(15)(16)K 二 0.003Gexp(k551(T+030),CX=0.(21)式中，C +和 C -為融土和凍土的容積熱容量；入+和入-為融土和凍土的導熱系數；Tp和Tb為劇烈相變區的上下溫度邊界； 衰減，其初始條件和邊界條件為：K為土體的導水系數，在劇烈相變區內隨溫度呈指數7(x. v.O) = 0(卞丁) T(

16、xt) = Fr(r). 丁(兀厶.r) =芮wox(17)(18)(19)(20)上述問題為溫度、水分影響下以穩態溫度場為初始溫度條件的瞬態熱傳導問題， 如果再加 上外力邊界條件和應力場的各種方程， 特別是用我們提出的本構模型作為物理方程， 則可 以進行溫度、水分影響下的應力計算， 即進行初步的水、熱、力三場的耦合計算。渠道靜力平衡方程為：(22)渠道靜力平衡方程為Zb = 0-幾何方程為&=，其中 d0 0 c0 dxdvc:0a5o 2 0丙dxe二00 0d c6二dy o物理方程(木構方程)為9HD倍(23)(24)(25)上述方程即組成了描述凍土渠道內含有相變的溫度場、水分場、應力

17、場的三場耦合方程組。3。1。2定解條件和各參數的取值在計算中取冰的體積含量 ci為0。05， 土的彈性模量 Es為12。0 MPa，冰的彈性模量取 為4910。0 MPa, 土和冰的泊松比均取為 0。3。整個平面上的初始溫度為 -0。3C。 上、 下邊界溫度變化為：| 片(O = -11.2siii(nH 3120.0)(”).F()二 10.9亡皿4恥(d).7 二計算模型以不良工程地質條件下渠道基土為例進行數學描述5。渠道垂直截一斷面為平面問題，認為土質均勻，水熱蒸發耗熱量及其它勢場忽略不計。整個渠道模型寬為3。0 m,高為2。0 m,頂端渠堤寬為1。0 m,渠底寬0。5 m ,渠底與下底

18、邊距離為 1。0 m。模 型的位移邊界條件為：下邊界Y方向限制移動、X方向自由， 左、右邊界X方向限制移動、Y方向自由，上邊界的X和Y方向均自由。應力邊界條件為：左右邊界剪應力均為0,右邊界上X方向應力為0。模型上邊界受相當于渠道襯切板重量和與襯切板垂直的凍 結力的均布壓力 P=1400N,方向垂直于邊界線，如圖3所示。 為了對比，圖4給出了文獻計算時的外載荷施加位置圖。各相關參數的取值為C+=3076。7, C-=2279。0 ( kJ m-3 C-1),入 + = 5。53,入-=6。58 ( kJ m-1 h-1 C-1),圖3本文計算時外載荷施加位置圖4 文獻 5的外載荷施加位置Ki=

19、1 o 8X 1O-o。7, K2=0。0036exp(0。551(T+0。3),K3=3 o 6X 10-3 (m h-i),L=334。7 (kJ kg-i),=-0 o 30, p T = -0o 75 b T (C ),0。0664 exp(0 。0551 ),0。75,0o 3058 0o 596( 0o 30), 0o 75 0o 30,0o 3058,0o 30ouT TW T TT? ? -3o 1o3 計算結果分析與討論采用三角形網格劃分, 共分成 12484 個三角形網格單元和 6463 個節點, 為了對整個渠 道凍結過程中的溫度、位移、應力、應變進行仔細的觀察, 我們取了

20、較長的時間來進行計 算，時間單位為小時，總共計算2000o 0個時間，時間步長取為5o 0o (1)溫度場分析。在渠道截面上， 不同凍結時間的0C等溫線分布計算結果如下圖 5所示，將其與圖6的 文獻計算的結果比較可以看出，在幾乎相同的凍結時間，0 C等溫線分布相似圖5本文計算得到的0C等溫線隨時間的變化18 分別為 50, 100, 200, 300, 500, 800, 1000, 1500 h 時的數值圖6文獻5 0 C等溫線隨時間的變化17 分別為 80, 320, 560, 800, 1040, 1280 和 1600 h 時的數值凍結2000 h后，上邊界的負溫達到最大值 -10o

21、12C(如圖7所示)，這與文獻 計算得到 的凍結1536 h后上邊界負溫最大值-10o 7C基本一致。從圖7中還可以看出，從渠坡到渠底的溫度變化快, 也即溫度梯度大, 而在離渠坡較遠的渠堤下, 溫度分布呈一組幾乎 平行的直線, 基本上不受渠坡邊界溫度的影響。 從溫度場的分布和變化來看, 二維邊界、 截面幾何形狀、 底水補給以及地中熱流對溫度分布均有顯著影響。 從以上比較來看, 溫度 場的計算結果是令人滿意的。 (2) 位移場分析。 大量的工程實踐表明, 在低凍結速率和有 充足水分補給的條件下, 產生較強的凍脹作用。 在渠道底部和渠坡下部, 溫度較高的地 下水起到降低凍結速率和補給凍結鋒面水分的

22、雙重作用, 在渠底和渠坡下部發生較強烈的凍 脹。 圖 8 是 2000 h 時的凍脹位移 ( x2 + y2 )分布, 可以看出, 最大凍脹位移發生在渠堤頂 端和渠坡上部, 渠底和渠坡下部位移量也比較大, 與同時段的觀測資料 (圖 9)對比,可見 渠道截面上的位移分布規律是基本一致的。圖7 凍結2000 h 時的等溫線分布 圖8 凍結2000 h 時的等位移線分布 圖9 寧夏大新渠試驗段 (南北走向 ) 實測表面位移分布(3) 應力場分析。 圖 1015 是渠道凍結 2000 h 后,各方向的應力和應變等值線分布圖。 在 渠底和渠坡下部的凍脹力都比較大, 在渠坡與渠底的交界點處出現了應力集中,

23、 最大凍脹 力也在此出現： 最大d x為3。943X 105 N/m2,最大r y為0。85X 105 N/m2,最大t xy為1。376 X 105 N/m2從圖中可以看出， 凍脹應力沿渠道基土深度方向的分布狀態。在渠道基土上部渠堤處， 由于外載荷圖10凍結2000 h的等d x線分布圖11凍結2000 h的等 x線分布圖12凍結2000 h的等dy 線分布的約束作用， 基土中的凍脹應力呈衰減趨勢。 隨著深度的增加， 外載引起的變 形逐漸減弱而由凍脹引起的變形逐漸增強，在超過某一深度后基土中的應力又會隨著深度的增加而增大， 而使基土應力隨深圖13凍結2000 h的等 y線分布圖14凍結200

24、0 h的等t xy線分布圖15凍結2000 h的等丫 xy線分布度的增加呈現衰減-增大-衰減的規律，這與文獻 5中計算得到的結論相同。從以上各圖還可以看出，應力和應變存在著一一對應關系， 等應力線和等應變線的變化趨勢相同說明我們所用的彈性計算模型是正確的。為將 計算結果看得更清楚一些， 我們選定了平面內的九個點作為研究對象， 分別來看各點的溫 度、位移、應力、應變隨時間的推移而產生的變化。圖 16 所選點的位置示意圖圖 17 各點的溫度隨時間變化曲線圖 18 各點 y 方向的位移隨時間變化曲線從圖 1620 可以看出， 在離渠道中心線 (右邊線 )越遠處， 溫度變化得越劇烈， 在靠近渠 坡處各

25、點的位移變化比較明顯。 另外， 由于左右兩邊界上位移條件的限制， 各個點水平 方向的位移不是很大， 而在垂直圖 19 各點 y 方向的應力隨時間變化曲線 圖 20 各點 y 方向的應變隨時間變化曲線 方向的位移變化則比較大。 靠近渠坡與渠底的交界處， 應力變化比較明顯， 量值也比距 渠堤較近處點的應力值大， 表明在渠坡和渠底的拐角處出現了應力集中，靠近這個區域點的應力變化比較劇烈。 而同一點處在同方向的應力和應變值， 在變化趨勢上是相同的， 它 們之間存在一一對應關系。3。 2 路基凍結三場耦合數值算例3。 2。1 控制微分方程控制微分方程同上面算例一樣， 連同方程 (22)(25)即組成了描

26、述凍土路基內含有相變 的溫度場、水分場、應力場的三場耦合方程組。3。 2。2 定解條件取路基橫斷面以二維問題來處理， 并取路基的一半作為計算區域。 設計算區域由 S1， S2, , S6六層均勻介質組成(圖21),其中水泥路面的半寬為 4。5m(包括路肩1。0 m), 路堤高 3。 0 m。 (1) 邊界條件。 根據青藏高原幾十年的觀測資料和附面層原理可得到固定 邊界上的邊界條件如下 17圖21 路基溫度場計算模型S1S6 為路基計算區域； B1B6 為路基區域邊界( )? n n ? = + +? U U0 0 01 2 3( ) sin 2 3 ， 8760 55f t T R t A t

27、x y B B B(27)? s T G ，n( ) 5 x， y B ， (28)( ) U 4 6 x， y B B ， (29)式中 T0 為下附面層底的初始年平均地溫 (地表以下 0。5 m 處的實測溫度值 )， 其取值 與文獻 17 的取值一致 (見表 2)， R0 為氣候變暖引起的下附面層底地溫增溫率， A0 為下 附面層底地溫振幅， n 為邊界界面的外法線方向， Gs 為計算區域下邊界地溫梯度， Gs=0。 03C/m。路基模型的位移邊界條為：下邊界B5限制y方向移動，x方向自由； 左邊界B4和右邊界B6均限制x方向移動，y方向自由；上邊界B1，B2和B3的x方向和y方 向均自由

28、。應力邊界條件為：左右邊界B4和B6上的剪應力均為0,右邊界B5上x方向應力為 0， 上邊界 B1、B2 和 B3 上無載荷作用， 因而在溫度變化過程中產生的應力全 部為熱應力， 土的熱膨脹系數取為 1。0X 105。表2 典型表面附面層參數值類型 天然地表 砂土邊坡 水泥路面附面層總溫度增量(C) 2。5 3。7 4。7(2) 初始條件。 根據邊界條件 (27)(29)和計算區域土體的熱參數及特性參數求解 Laplace 方程： e T e T 0， x x y y入入(30) 得到季節融化開始時計算區域溫度場的估計值。3。 2。3 參數的選取根據工程實際， 本算例的計算區域由六種不同的土質

29、組成(見圖 21 所示 )， 其熱參數及土壤特性參數的取值如下表 3 所示。 表中 C 為路基土的容積熱容 表3 計算用土體熱參數及特性參數土類參數S1砂礫與碎石土 S2草炭亞黏土 S3含泥炭亞黏土 S4含碎石黏土 S5含碎石亞黏土 S6弱風化基巖入u(J m-1 C-1)6875。1 3376。7 5110。8 2052。6 3124。2 7164。8入f(J m“ C-1)9396。0 4684。1 6249。1 3162。1 4376。2 9689。2Cu (J m-3 C-1) 2。183X106 2。258X106 2。342X106 2。166X106 2。446X106 2。78

30、4X106Cf (J m-3 C-1) 1o 694X106 1。781X106 1。932X06 1。539X106 1。932X06 2。192X1C6P d(kg m-3) 1800 1200 1400 700 1300 2250W(wt%) 10 24 22 56 24 5。 6Wu(wt%) 3 7 7 20 7 0。 6D 9。 35X10-6 1。 65X10-6 4。 66X10-6 1。 21X10-6 3。 73X10-6 3。 44X10-6量，入為路基土的導熱系數，下標f和U分別表示路基土體的凍、融兩相狀態，pd 為土體的干容重 (干密度 )， W 和 Wu 分別為路基

31、土的重量含水量和重量未凍水含量，D為路基土的水分擴散系數， 一般認為土的 x 方向和 y 方向的水分擴散系數相等，即 Dx=Dy。未凍水含量與溫度的關系根據下式確定 5：( )Co 0644e。0551 , Co 75,0。3058 0。596 0。30 ， 0。75 0。30，0。 3058,0。 30。TwTT TT0? (31)3。 2。4 計算結果與分析計算采用三角形網格單元對整個路基區域進行網格劃分， 共分為 11800 個單元， 6061 個節點。 (1) 溫度場分析。 根據邊界條件 (27)(29)和計算區域土體的熱參數及特性參數求解 拉普拉斯方程 (30)， 得到季節融化開始時

32、計算區域溫度場的估計值，如圖 22 所示。在上述拉普拉斯方程解出的初始溫度分布基礎上，利用溫度條件表達式(27)中R0=0。0時作為上邊界的溫度條件逐時段求解方程 (10)， 直到年變化層以下溫度場基本保持穩定， 而且年變化層以上相同位置的溫度值在同一時間逐年相同時為止，此時， 各個節點的溫度值作為三場耦合計算時的初始條件， 溫度分布如圖 23 所示為了對整個路堤凍結過程中的溫度、應力、應變和位移進行詳盡的觀察，我們取了較長的時間來進行計算， 時間單位為小時 (h)， 總共計算 2000。0 個時間， 時間步長取為 1。 0。 圖 2427 是凍結不同時間后路基的溫度等值線變化圖。可以看到，由

33、于路基上邊界有隨時間變化的溫度邊界條件， 下邊界有熱流交換， 所以在凍結初始段路基內部的溫度變化 比較劇烈，最低溫度在凍結 50 h后就由初始時刻的-9。229C下降到-9。734C,在500h達 到了-10。24C；最高溫度也由初始時的0C下降到了 500 h的-3。988C。由于上邊界是周期變化的溫度條件， 所以經過一段時間后路基平面內的溫度又有所回升，最低溫度由500 h 的-10。24C上升到了 1000 h 的-9。 821 C和 2000 h 時的-9。 145C。 此外， 從整 個溫度場的變化可以看出， 溫度變化比較劇烈的區域都是在路基上部距離天然地表約5。 0m6。 0 m 的

34、深度范圍內， 在此深度以下的大部分區域里溫度變化都不明顯。在地表下約1。 0 m3。 0 m 的深度范圍內存在一個溫度變化更為頻繁、劇烈的區域，這就是通常所說的劇烈相變區域， 在此區域內水、 冰之間由于凍融而產生相變轉化。 這些現象與實際工程 應用中的觀察結果一致，一般受氣溫變化影響的地表深度也是在5。0 m10。0 m以內，比這個深度更深的區域基本不受地表氣溫變化的影響。(2) 位移場分析。 由于本節對路基進行耦合計算時模型不受外載荷的作用，所以在計算中所產生的應力、應變、位移全部由凍脹和融沉產生。從上面的圖 2829 中不同時刻的位移變化情況看， 整個凍結過程中產生的水平方向 (x 方向

35、) 的最大位移約為 0。0002853， 垂直方向 (y 方向 ) 的最大位移約為 -0。000472， 這說明凍脹產生的水平方向變形稍小于垂 直方向的凍脹變形。 (3) 應力場分析。 圖 3031 是出現最大應力時的云圖， 整個凍結過程 中產生的水平方向(x方向)的最大應力約為1。81 x 105 Pa垂直方向(y方向)的最大應力約 為0。4974X 105 水平方向的應力是垂直方向應力的三倍以上， 原因與在計算時限制了水平方向的位移有關。從圖中還可以看出， 水平方向的最大應力與垂直方向的最大應力是在同一時間段出現的，這說明由于隨著邊界溫度條件的周期性變化， 凍土路基內的應力、 變形也隨之發

36、生相應的變化， 并且這種外溫度 邊界條件的影響區域也是距離天然地表約6。0 m深度以內，在深于6。0 m的區域則影響較小。圖 3235 分別是凍結 2000 h 時各應變的云圖。 從計算結果可以看出， 距離路基中 心線(右邊界)越近、越是靠近路堤和路坡處，無論是溫度變化， 還是應力、位移、應變的變化， 都比路基內其它區域要劇烈一些， 這正是因為這些區域更容易受地面溫度變化的影 響所致。 在整個路基區域的不同凍結時刻， 應力等值線與相應方向的應變等值線圖形非常 相似， 說明應力與應變還是一種一一對應關系， 原因還是由于我們所使用的本構模型是一 種彈性模型。 我們對凍土路基的三場耦合計算結果與前人

37、的計算結果以及工程實際中的觀 察結果相符合， 規律相似， 說明我們所做的三場耦合分析是可靠和可信的。此外，我們 從路堤頂端開始在距離路堤中心線圖22 拉普拉斯方程解出的初始溫度分布 圖23 三場耦合計算時的初始溫度分布 圖24 凍結50 h 時的等溫線分布 圖25 凍結500 h 時的等溫線分布 圖26 凍結1000 h 時的等溫線分布 圖27 凍結2000 h 時的等溫線分布 中國科學： 物理學 力學 天文學 2010 年 第 40 卷 第 6 期 769圖28凍結1000 h時x方向的等位移線分布圖29凍結1000 h時y方向的等位移線分布圖30凍結500 h的等x線分布圖31凍結500

38、h的等y 線分布圖32凍結2000 h的等& x線分布圖33凍結2000 h的等y 線分布朱志武等： 基于損傷的凍土本構模型及水、熱、力三場耦合數值模擬研究770圖34凍結2000 h的等t xy線分布圖35凍結2000 h的等丫 xy線分布（右邊界 ）3。0 m 處、沿平面內每隔 1。5 m 深度選定八個點作為研究對象， 分別看它 們的溫度、位移、應力、應變隨時間的推移而產生的變化。 而且由于整個凍融過程在天然 地表以下 10。0 m 深度后各物理量的變化比較小， 所以我們所選擇的點都在距離路堤頂端 15。0 m 以上。 所選定 8 個點的位置如圖 36 所示。圖 3740 即是以上各個點各

39、量值的計 算結果比較： 可以看出， 修筑路基后由于路基內部的儲熱作用， 使得凍土路基垂直方向的 最大融化深度在竣工的前幾年迅速增大， 以后逐漸趨于穩定。 在五、 六年之后凍土路基溫 度場的影響已經較小， 并且在較深的路基內部， 溫度場的變化越來越緩慢， 這與實際工 程實踐中的觀察結果相吻合。圖 36 8 個特殊點的位置示意圖圖 37 溫度隨時間變化曲線圖 38 y 方向的位移隨時間變化曲線圖 39 y 方向的應力隨時間變化曲線圖40 y 方向的應變隨時間變化曲線在距離路基頂端越近處點的溫度、位移、應力、應變的變化都較劇烈，說明隨著深度的不同， 外界條件的化對淺層地表的影響比較大。 這也從另外一

40、個角度說明， 在實際施 工過程中， 如果把路堤修得越高， 天然凍土上限越淺， 凍土路基中的最大融化深度也越 淺， 越有利于保護天然凍土， 使其受外界施工的影響越小。 但是在實際的施工過程中采 用單純提高路堤高度的方法也是很不經濟的， 另外一種辦法就是對路堤上表面的覆蓋材料 進行選擇， 選擇那些傳熱性能差的材料鋪設路基， 也可以在路堤上鋪設保溫層來防止外界 溫度變化對路基內部天然凍土的影響， 以保證凍土的天然上限穩定。4 結論本文針對土體凍融過程中水、熱、力三場耦合這一重要的難點問題， 用細觀力學的分 析方法， 將凍土內部各組分的力學行為與宏觀的凍土力學行為結合起來進行研究， 在宏觀 和細觀之間

41、建立了一種橋梁，在研究方法上為深化凍土力學的研究提出了一種新的分析思 路。 (1) 建立了處于某一特定溫度情況時土、冰的彈性常數與整個凍土體彈性常數之間的關 系， 并在此基礎上， 進一步將損傷力學的方法應用于凍土力學研究中， 建立了考慮土體 損傷的凍土本構模型。 (2) 渠道凍結三場耦合數值計算結果與前人計算結果及其實測結果相 符合， 規律相似。 計算結果表明二維邊界、 截面幾何形狀、 底水補給以及地中熱流對溫度 場分布均有顯著影響； 隨著深度的增加， 外載引起的變形逐漸減弱而由凍脹引起的變形逐 漸增強，在超過某一深度后基土中的應力又會隨著深度增加而增大。 (3) 凍土路基的三場耦 合計算結果

42、表明， 在整個路基上部約距離天然地表 5。06。0 m 的深度范圍內溫度變化比 較劇烈， 并在地表下約 1。03。0 m 的深度范圍內存在一個溫度化更為頻繁、劇烈的相變 區， 在此區域內水、冰之間由于凍融而產生相變轉化。 在距離路基中心線 (右邊界 )越近、 越是靠近路堤和路坡處， 因為更容易受地面溫度變化的影響所致， 所以無論是溫度變化， 還 是應力、位移、應變的變化， 都比路基內其它區域要劇烈一些。參考文獻Fremond M, Mikkola M. Thermomechanical of freezing soil. In ： Proceedings of the Sixth Intern

43、ational Symposium on Ground Freezing.Rottendam : A. A. Balkema, 1991.57 -66Gary P, Jason W. A model for freeze-thaw cycling of the South Pole Station exterior sheathing. Cold Regions Sci Technol, 2008, 52 :254 -61.Konrad J, Duquennoi C. A model for water transport and ice lensingin freezing soils. W

44、ater Resour Res, 1993, 29 :3109 -3124Dennis E, Pafahl A. Frost Action. Section 3 of Roads and Airfields in Cold Regions. Edited by Teds Vinson. Washington ： Washington State Uni versity Press, 1996. 51 -85An W D, Wu Z W, Ma W, et al. Interaction among Temperature Moisture and Stress Fields in Frozen Soil (in Chinese). Lanzhou ： Lanzhou University Press, 1989. 78 -95Foriero A, Ladaryi B, Dallimore R, et al. Modelling of deep seatedhill slope creep in permafrost. Can Geotech J. 1998, 35 :560 -578Su B, Li N, Quan X J. The numerical