地下水流影響下含鹽飽和粉砂地層中聯絡通道凍結溫度場的特征分析

李 飛, 王長虹, 張海東, 董濟涵

(上海大學力學與工程科學學院, 上海 200444)

人工地面凍結法是將土體周圍地層凍結, 加固土體、隔絕地下水, 并把對鄰近構筑物的影響降低到最小程度的一種工法, 已被廣泛應用于地鐵工程[1-6].地鐵聯絡通道施工是地鐵建設的薄弱環節, 常處于復雜的地質環境, 如在沿海地區修建地鐵聯絡通道, 存在水流和鹽分等不利地質條件, 既要考慮聯絡通道本身的安全, 又要考慮對周邊建筑物的影響, 若處理不當勢必引起嚴重的工程事故.

關于凍結過程中的多場耦合分析已經有了很多研究成果.在理論和數值模擬方面,Harlan[7]基于非飽和土中水分遷移與凍結土層中未凍水的水分遷移理論, 提出了水-熱耦合概念.Lai 等[8]研究了飽和粉土水-熱-力相互作用過程的理論模型.Grenier 等[9]研究了凍融系統的地下水流動和熱交換.岳漢森[10]在Harlan 模型的基礎上, 考慮土體凍融過程中鹽分的對流和擴散作用, 建立了水-熱-鹽耦合運移數學模型.賴遠明等[11]根據傳熱學、滲流理論及凍土力學提出了考慮相變的溫度場、滲流場和應力場耦合問題的數學模型及其控制方程, 并應用于巖土工程.肖澤岸等[12]建立了含氯化鈉鹽土層凍結過程中水-鹽遷移及變形的計算模型, 該模型反映了含鹽土體層在凍結過程中的溫度、水分、鹽分及變形規律.Xu 等[13]考慮溶質的擴散與對流, 基于Darcy 滲流定律、能量守恒定律和Fick 擴散定律提出了鹽漬土的水-熱-鹽耦合方程.

在室內試驗方面, Pimentel 等[14]在滲流條件下開展了人工凍結土的大型室內試驗研究.邴慧等[15]通過對青藏鐵路沿線不同土質在不同含鹽量和含水量下的凍結溫度試驗, 得出土層凍結溫度隨含鹽量的增加而降低, 隨含水量的增加而升高的結論.孫立強等[16]通過室內試驗得出了土體熱物理參數隨溫度場的變化規律.

在現場監測方面, 林萍等[17]以上海長江隧道的聯絡通道為工程背景, 提出了一種凍土壓力的現場監測方法.蘇文德[18]針對國內首條海底地鐵隧道的凍結法展開研究, 揭示了海底隧道聯絡通道凍結溫度場的變化規律.Yan 等[19]采用現場監測、理論解析和數值模擬相結合的方法, 分析了廣州地鐵3 號線聯絡通道的凍結壁厚度與平均溫度.郜新軍等[20]基于現場試驗,對飽和粉質黏土地層聯絡通道凍結法施工中的溫度場變化進行了研究, 對溫度、凍結壁、凍脹和地表沉降的發展規律展開了分析.

上述研究從基礎理論、數值計算、室內試驗和現場監測等方面對凍結法進行了總結, 對凍結法的溫度場、滲流場和鹽分場隨時間的變化進行了深入分析.但類似南通地鐵聯絡通道凍結法在地下水流-鹽分共同作用下的機理尚不明確.為此, 本工作以南通城市軌道交通1 號線永興大道站—深南路站區間4#聯絡通道工程為背景展開研究.首先, 建立了水-熱-鹽三場耦合的物理數學方程; 然后, 利用有限元軟件COMSOL Multiphysics 模擬含鹽土體的凍結模型試驗,并分析對比試驗結果; 最后, 基于背景工程的特點, 開發了多物理場數據監測設備, 對比分析監測數據與三場耦合模擬數據.

1 控制方程

為了研究含鹽土層在人工凍結過程中的水-鹽遷移和水-冰相變過程, 考慮水流影響, 建立水-熱-鹽三場耦合方程.引入以下假設: ①土體骨架和冰晶體不可壓縮; ②土體為各向同性的均勻介質; ③水分在遷移的過程中符合Darcy 定律.

1.1 相變溫度與未凍水含量

鹽分的存在降低了含鹽地下水的凍結溫度, 進而影響土層的凍結過程.通過Clausius-Clapeyron 方程可得到相變溫度的降低幅度.通過適當簡化, 凍結溫度Tf的變化[21]為

式中:Tf為土體凍結溫度(K);T0為純水凍結溫度(K);R是氣體常數,為8.314(J·mol-1·K-1);ρw為水的密度(kg/m3);L為水-冰相變時釋放的潛熱值(kJ/kg);c為鹽分濃度(mol/L).含鹽土體層的凍結溫度Tf也可以表示為

在孔隙溶液中, 如果有鹽分的存在, 不僅會影響土體凍結溫度, 還會影響溶液的成分.在不含鹽分的土體層中, 孔隙中未凍水的飽和度Sw0[9]為

式中:Sr為殘余含水量;W為與土性有關的參數.

當土體層含有鹽分時, 含鹽孔隙溶液的飽和度Sl[12]為

式中:A1,A2和A3是與土性有關的系數;M,m是與孔隙有關的系數.

孔隙中未凍水的飽和度為

式中:V?為鹽分的表觀摩爾體積.

1.2 水分場控制方程

水分的變化符合Darcy 定律, 根據質量守恒定律, 孔隙溶液中水分的質量變化[9]表示為

式中:n為孔隙率;Si,Sw分別為冰和水的飽和度,Si+Sl= 1;U是水流速度,U=Kw為滲透系數.隨著溫度的降低, 冰晶將堵塞孔隙, 導致土體的滲透性降低, 則孔隙溶液中土體的滲透系數[22]為

式中: flc2hs 為平滑的Heaviside 函數;δT為相變區間;kr為殘余滲透系數.不考慮冰晶的壓縮性, 式(6) 可以表示為

考慮水的壓縮性, 引入水的壓縮系數[9]

則式(8) 可以表示為

當溶液為稀溶液時, 濃度對未凍水含量的影響很小,Sl可以表示為未凍水的含量, 則水分場的控制方程可以簡化為

1.3 溫度場控制方程

在傳熱過程中, 熱傳遞包括3 個方面, 分別為熱傳導、熱對流和熱輻射.由于地下水的流動, 對流所引起的熱傳遞不可忽略, 地下水流速度越快, 對傳熱的影響就會越大.熱輻射的量很小, 可以忽略不計.此外, 水-冰相變將釋放大量的潛熱, 能量守恒方程除了熱傳導和熱對流, 還包括相變項.溫度場的能量平衡方程[9]為

式中:ρs、ρl和ρi分別為土體骨架、孔隙溶液和冰晶的密度, kg/m3;Cs、Cl和Ci分別為土體骨架、孔隙溶液和冰晶的比熱容, J/(m3·K);λs、λl和λi分為別為土體骨架、孔隙溶液和冰晶的熱傳導系數, W/(m·K).

1.4 鹽分場控制方程

土體中溶質的運移包括3 種方式, 分別為對流、分子擴散和機械彌散[23].在含鹽地下水的流動過程中, 主要涉及對流和分子擴散.鹽分的變化主要來源于兩個方面, 一方面是對流引起的鹽分遷移, 對流指的是鹽分溶解于液態水中而遷移的過程, 在凍結過程中, 主要指的是未凍水在流動過程中引起鹽分的變化; 另一方面是分子擴散引起的鹽分遷移, 分子擴散指的是溶液中溶質分子的熱運動, 可在濃度梯度的條件下發生, 不受溶質分子運動的制約.根據質量守恒定律, 鹽分場的遷移方程[12]為

式中:Ms為溶質的相對分子質量;D為溶質的擴散系數(m2/s),D=D0Slξ,D0為溶質在未凍區的擴散系數.

2 多場耦合分析驗證

參考經典的鹽漬土凍結模型試驗, 下面采用數值計算方法模擬凍結模型試驗, 對比分析模擬結果與試驗數據, 驗證多場耦合分析方法的有效性.

2.1 模型試驗

應賽等[24]在鹽漬土凍結過程溫度變化的研究中, 對土樣在不同濃度下進行了室內凍結模型試驗.如圖1 所示, 將高35 mm, 直徑35 mm 的土樣放入金屬容器中, 采用低溫冷浴控制裝置對土樣進行降溫凍結.將溫度探頭插入金屬容器內的土樣, 冷浴溫度由18?C 降到-15?C,土樣的降溫速率為(0.3±0.05)?C/min, 頂板為絕熱材料.

使用COMSOL Multiphysics 有限元軟件自定義偏微分方程模塊求解土樣在不同鹽分濃度下, 溫度隨時間的變化規律.在多場耦合分析中, 土樣的鹽分濃度分別取0、0.076、0.305 和0.609 mol/L.通過數值計算驗證土樣凍結溫度和凍結過程隨鹽分的變化規律.

2.2 模擬與試驗結果對比分析

土樣在不同鹽分濃度條件下的模型試驗結果如圖2 所示.鹽分濃度越高, 土樣降溫凍結至-15?C 所需的時間越長.凍結過程經歷4 個階段: 第一階段為過冷階段, 孔隙溶液沒有開始轉化成冰, 此時的最低溫度稱為過冷溫度, 也是冰開始產生的溫度; 第二階段為上升階段, 大量冰開始產生, 釋放大量潛熱, 土樣溫度略微回升; 第三階段為穩定階段, 冰水相變所釋放的潛熱與外界交換熱量相平衡, 溫度不變; 第四階段為下降階段, 當土體中的自由水大部分形成冰后,將不會產生潛熱, 因此土體溫度快速下降.在凍結前期, 不同鹽分的土樣之間溫度變化不明顯,這是因為凍結前期并沒有形成冰, 鹽分主要對土體凍結溫度產生影響.數值模擬結果和模型試驗數據契合度較高, 誤差在0.5?C 之內.

圖2 不同鹽分濃度下溫度隨時間的變化Fig.2 Temperature variation over time at different salt concentrations

利用有限元軟件COMSOL Multiphysics 多場耦合分析方法能夠有效模擬含鹽土體的凍結過程.已有研究[9,25]中, 含有水頭的人工凍結法的模擬已經較為成熟, 因此可以將水-熱-鹽三場耦合分析方法應用于水流-鹽分共同作用下的聯絡通道的凍結施工模擬.

3 工程應用

以南通城市軌道交通1 號線永興大道站—深南路站區間4#聯絡通道為工程背景, 開發了多物理場數據監測設備.根據聯絡通道的實際尺寸, 建立聯絡通道凍結施工二維數值模型進行模擬分析.對比監測數據與模擬數據, 分析鹽分與水頭對凍結壁平均溫度和凍脹變形的影響.

3.1 工程背景

南通城市軌道交通1 號線永興大道站—深南路站區間起于永興大道與長泰路交口的永興大道站, 呈直線形和曲線形下穿長泰路和永和路, 進入深南路站, 此區間共包含4 處聯絡通道.如圖3 所示, 4#聯絡通道及泵房位于永和路下方, 與永興佳園60 幢距離最近處約31 m, 其中心里程為左線SK13+610.373, 隧道中心標高-18.572 m; 右線XK13+631.906, 隧道中心標高-18.584 m.聯絡通道線間距14.210 m, 聯絡通道埋深21.70 m, 聯絡通道所在土層為飽和粉砂層.聯絡通道采用人工凍結法施工, 凍結壁厚度的設計要求為不小于2.2 m, 平均溫度低于-12?C.

3.2 多物理場數據監測設備

傳統人工凍結施工中僅布置溫度監測孔, 且傳感器位于凍結管內, 并未接觸土層, 無法實現同一位置多物理場數據的監測.因此, 開發了一體化多物理場數據監測設備如圖4(a) 所示,監測管一端置于凍結壁內部, 另一端連接監測通訊設備.監測設備是一種防水集成監測管裝置, 包括端部鏤空鋼管、普通鋼管、隔水部件和傳感器.隔水部件的兩端分別通過套筒螺紋連接端部鋼管和普通鋼管, 傳感器設置于端部鋼管內, 用于采集土層多物理場數據.3 種傳感器分別為BSIL-T2 型溫度計、TDKYJ30 型孔壓力傳感器、TR-HTS03 三探針溫度-濕度-鹽分一體化傳感器.3 種傳感器用于監測凍結土體中的孔隙水壓力、溫度、鹽分和含水量等數據.傳感器測得的數據通過DataTaker DT800 數據采集器進行記錄和傳輸(見圖4(b)), 為凍結法施工提供多物理場監測數據.

圖4 多物理場數據監測設備和數據采集器Fig.4 Multi-physical field data monitoring equipment and data collector

3.3 凍結數值模型

針對南通市軌道交通1 號線永興大道站—深南路站區間4#聯絡通道建立數值計算模型.凍結管的布置如圖5 所示, 以凍結管的對稱軸為中心, 橫坐標為與對稱中心的距離, 縱坐標為距離地面的深度.M1～M4為工程監測點, 僅能測得溫度;M5為一體化多物理場數據監測點, 能夠測得溫度、水壓力、鹽分和含水率等數據.人工凍結管對稱布置在聯絡通道的開挖外輪廓, 通過分析監測點溫度數據, 發現兩側溫度監測點的數值存在差異, 因此判斷存在地下微水流, 方向為右向左流動.然而, 除了地下水流的作用, 鹽分和導熱系數也會影響凍結時間.根據溫度測試數據, 聯絡通道所處土層的初始溫度約為23.7?C, 凍結管的側壁是冷凍系統源, 通過鹽水循環提供冷區溫度, 鹽水溫度的變化如圖6 所示.通過對凍結管兩側溫度監測點數值的試算, 設置滲流水頭h分別為0.02、0.05 和0.08 m, 導熱系數λs分別為1.7、1.9 和2.1 W/(m·K).飽和粉砂層的多物理場力學參數如表1 所示.

表1 多物理場計算參數Table 1 Multi-physical field calculation parameters

圖5 凍結管和監測點布置Fig.5 Layout of freezing pipes and monitoring points

圖6 鹽水溫度變化曲線Fig.6 Temperature change curve of brine

3.4 結果分析

下面通過對模擬結果與監測數據進行對比分析, 討論水頭、鹽分和導熱系數對凍結過程的影響.基于模擬結果分析凍結壁的平均溫度和厚度.

3.4.1 模擬結果與監測數據對比

隨著凍結管溫度的降低, 土層溫度也隨之發生變化.在相同環境下, 各點溫度的變化趨勢一致.圖7(a) 分別為不同水頭和不同導熱系數下,M1點溫度的變化.由圖可知, 在凍結前期,土層的溫度變化較快, 因為凍結管的溫度與土層溫度相差較大.對比導熱系數不變的幾條曲線可知, 隨著水頭的增大, 溫度降低越來越慢, 這是因為在地下水流動過程中, 水流帶走了一部分能量.當水頭保持不變時, 隨著導熱系數的增大, 溫度降低越來越快.當導熱系數不變時, 隨著水頭的增大, 溫度降低越來越快(見圖7(b)).隨著水流從右邊流到左邊, 左邊的監測管獲得了更多熱量, 并且當凍結壁形成后, 左邊受水流的影響會更小, 因此溫度下降更快.水頭使凍結壁整個向左遷移, 這對凍結時間有較大影響.分別對比圖7(a) 和(b) 中模擬數據與監測數據,得出當h=0.05 m,λs=1.9 W/(m·K) 時, 可以獲得擬合效果相對較好的計算結果.

圖7 不同水頭和導熱系數下溫度的變化曲線Fig.7 Temperature variation curves under different water heads and thermal conductivities

圖8 為M3～M5點監測與模擬數據溫度的對比曲線.經過55 d 的凍結,M3和M4點溫度變化較M5點相對平緩, 而M3點溫度高于M4點, 這是因為水流從右往左流, 帶走了一部分熱量.模擬結果與監測數據相似度較高, 模型符合實際情況.M5點水壓力的變化如圖9 所示, 在凍結前期, 水壓力基本保持在220 kPa 左右; 大約凍結400 h 時, 水壓力發生了跳躍, 達到240 kPa, 從水分場的平衡方程可以得出, 在水變成冰的過程中, 短時間內流進該點的水分不會有很大變化, 滲透系數會急劇減小, 因此水壓力突然增大; 隨著繼續凍結, 水壓力慢慢趨于穩定, 最后保持在227 kPa 左右, 這是因為隨著不斷的凍結, 外部流進該點的水分會不斷減少,水流穿過孔隙慢慢趨于穩定, 即水壓力也會恢復穩定.

圖8 M3、M4 和M5 點監測與模擬數據溫度對比Fig.8 Temperature comparison of monitored and simulated data at points M3, M4 and M5

圖9 M5 點水壓力監測與模擬值對比Fig.9 Comparison between monitored and simulated water pressure at M5 point

M5點鹽分濃度的變化如圖10 所示.在凍結前期, 濃度幾乎不變, 保持在0.1 mol/L; 大約凍結300 h 時, 濃度開始下降, 從0.1 mol/L 降到0.08 mol/L, 這是因為隨著該點附近被凍結,水流向凍結處遷移, 導致鹽分也隨之遷移, 此時鹽分濃度會有下降的趨勢; 隨著繼續凍結, 由于鹽分數據采集器的原因, 凍結后并未測出監測值, 而模擬值在凍結380 h 左右, 濃度發生較大跳躍, 從0.085 mol/L 增大到0.17 mol/L, 在水-冰發生相變的過程中, 水分慢慢變成冰晶, 未凍水減少, 濃度突然增大, 濃度的增大會反過來抑制未凍水繼續被凍結; 隨著凍結繼續, 濃度慢慢趨于穩定, 最后大致保持在0.15 mol/L, 當水壓力趨于穩定時, 未凍水的含量趨于穩定, 濃度也將保持穩定.

圖10 M5 點鹽分濃度監測與模擬值對比Fig.10 Comparison between monitored and simulated salty concentrations at M5 point

3.4.2 鹽分和水頭的影響

圖11 為M1點有、無鹽分和水頭差的溫度變化.由圖可知, 含鹽且有水頭的溫度變化最慢, 而不含鹽且無水頭的溫度變化最快.當凍結1 320 h 時, 含鹽有水頭工況下降到-10?C;不含鹽無水頭工況達此溫度只需1 080 h; 含鹽無水頭工況達此溫度需1 175 h, 比不含鹽無水頭工況延遲了95 h; 不含鹽有水頭工況達此溫度只需1 255 h, 比不含鹽無水頭工況延遲了175 h.從上述數據可以得出, 鹽分和水頭能夠延遲該點到達某一溫度的時間, 進而影響總體的凍結時間.根據溫度場的控制方程, 能夠看出水頭通過對流的方式帶走了一部分能量, 導致溫度下降更慢, 而鹽分主要影響液體的凍結溫度和熱物理參數, 導致溫度降低更慢.

圖11 水頭和鹽分作用下M1 點的溫度變化Fig.11 Temperature change under water head and salt of M1 point

3.4.3 凍結壁厚度和平均溫度

判斷聯絡通道是否達到開挖的條件主要依據凍結壁的厚度以及凍結壁的平均溫度, 圖12為凍結不同時間測線1 和測線2 的溫度分布.設計要求凍結壁的厚度為2.2 m, 虛線為設計要求范圍, 通過積分法求該范圍內的平均溫度為凍結壁的平均溫度.由圖可知, 溫度的分布形式為“W” 形, 右側溫度整體略高于左側溫度.

圖12 不同凍結時間溫度分布Fig.12 Temperature distribution with different freezing time

圖13 為上側、下側、左側和右側凍結壁平均溫度和厚度隨時間的變化.表2 為各側凍結壁溫度達到-12?C 和厚度達到2.2 m 所需時間.各側凍結壁平均溫度達到-12?C 所需時間分別為32、42、45 和50 d, 厚度達到2.2 m 所需時間分別為41、44、44 和48 d.因此, 通過數值模擬可得, 滿足設計要求的最少凍結時間為50 d.

表2 凍結壁厚度達到2.2 m 和平均溫度達到-12 ?C 所需時間Table 2 Time required for the thickness of the freezing wall reaching to 2.2 m and the average temperature of -12 ?C d

圖13 凍結壁厚度和平均溫度隨凍結時間的變化Fig.13 Variation of freezing wall thickness and average temperature with freezing time

由圖12、圖13 和表2 可知, 凍結壁上壁的厚度和平均溫度發展最快, 凍結41 d 就達到了設計要求; 凍結壁右壁的厚度和平均溫度發展最慢, 凍結50 d 才達到設計要求; 凍結壁的厚度和平均溫度與水頭密切相關, 靠近水頭一側的凍結壁厚度比另一側要薄, 平均溫度要高; 凍結壁厚度在凍結15～25 d 時發展速度最快, 此時正處于凍結壁交圈的時間段.

4 結 論

本工作基于多孔介質的熱傳導、滲流和鹽分的遷移理論, 建立了耦合的物理數學模型, 利用有限元軟件COMSOL Multiphysics 模擬南通地鐵飽和粉砂地層中聯絡通道凍結法的施工過程, 得出了如下主要結論.

(1) 考慮鹽分對土層凍結溫度和未凍水飽和度的影響, 建立了含鹽土層的水-熱-鹽三場耦合方程, 使用COMSOL Multiphysics 數值分析軟件建立計算模型, 并在不同水頭和導熱系數條件下對溫度進行確定性分析, 對比監測和模擬數據, 驗證了本方法能較好反映含鹽土體的凍結過程.

(2) 在土體發生水-冰相變的過程中, 水壓力會發生突變, 隨著繼續凍結, 水壓力會趨于穩定.因為當水-冰相變時, 冰會快速填充孔隙, 導致水流很難通過孔隙, 造成水壓力增大, 凍結完成后, 水流減小, 趨于穩定, 壓力隨之保持穩定.

(3) 含鹽土層比不含鹽土層降溫速度更慢, 并且隨著鹽分濃度的不斷增加, 對土層溫度的影響增大, 達到同樣溫度需要的時間更長, 這是因為隨著鹽分濃度越來越高, 凍結溫度會不斷降低, 導致溫度的變化過程后移, 因此在規定的凍結時間內, 濃度越高, 溫度越高.鹽分濃度還能通過影響凍結溫度來影響未凍水的飽和度, 進而影響凍結壁的厚度.水頭主要對上游一側產生影響, 上游一側的凍結壁厚度小于下游一側凍結壁, 這是因為當水流動時, 帶走了一部分能量, 使得上游一側更難以凍結, 造成右側凍結壁厚度更薄.

(4) 測線部分的溫度呈“W” 形分布.右側凍結壁發展最為緩慢, 厚度達到2.2 m 需凍結48 d.凍結壁達到2.2 m 的同時平均溫度達到-12?C 需凍結50 d.