王 彬，榮傳新，程 樺

(1.安徽理工大學 土木建筑學院，安徽 淮南 232001；2.安徽理工大學 礦山地下工程教育部工程研究中心，安徽 淮南 232001；3.安徽理工大學 安全科學與工程博士后科研流動站，安徽 淮南 232001；4.中煤礦山建設集團有限責任公司 博士后科研工作站，安徽 合肥 230091)

人工地層凍結法由Friedrich Poetsch首次提出，最初運用于富水軟土地層的立井掘砌施工過程中。在該工法的實施過程中，外部制冷機組首先將低溫冷媒降低至負溫（鹽水或者酒精的溫度一般被降至–40℃，液氮的溫度可以降至–190℃）；隨后，低溫冷媒被泵送至每根凍結管，在凍結管中的低溫冷媒通過管壁與周圍的土層進行熱量交換，以實現(xiàn)對土層的凍結。由于該工法對周圍環(huán)境的影響較小，其形成的凍結壁具有良好的封水性，作為臨時支護結構具有較高的強度，該工法已經(jīng)逐漸發(fā)展成為復雜地層地下工程施工的主要工法之一。近些年，在中國的天津、廣州、深圳、寧波等沿海城市的地鐵隧道建設中，凍結法廣泛應用于地鐵聯(lián)絡通道施工、盾構進出洞口土層加固等工程。受海相沉積環(huán)境、海水浸漬和潮汐影響，沿海城市的地層具有地下水礦化度高、流速大等特點，采用凍結法施工時，凍結溫度場的分布規(guī)律與無滲流場作用時存在較大的差異。

凍結壁的厚度及平均溫度等重要參數(shù)與凍結溫度場的分布規(guī)律存在著密切的聯(lián)系，因此，對凍結溫度場計算理論的研究尤為重要。人工地層凍結是一個隨時間變化的過程，其過程描述需要采用瞬態(tài)導熱理論。但由于數(shù)學計算上的困難，人工地層凍結瞬態(tài)溫度場解析解的求解至今主要針對單管凍結。而考慮滲流場對凍結溫度場的作用后，凍結鋒面的形狀、大小隨著滲流速度、凍結時間、介質的性質而發(fā)生變化，并且凍結鋒面作為相變邊界需要考慮冰水相變的過程，因此求解瞬態(tài)解析解的難度很大，目前僅Victor進行了單管凍結溫度場的相關推導，且是將凍結鋒面的形狀簡化成圓形，其計算結果與實際凍結溫度場存在一定差距。

人工地層凍結在進入穩(wěn)定凍結階段后，熱傳導速率較慢，其溫度場與穩(wěn)態(tài)溫度場非常接近。因此，可以在該階段假定熱傳導達到穩(wěn)定狀態(tài)，并采用穩(wěn)態(tài)溫度場近似等效瞬態(tài)溫度場。當?shù)貙又写嬖跐B流場作用時，進入穩(wěn)定階段后，在滲流速度及水溫不變的條件下，凍結鋒面處水流與被凍土體之間的對流傳熱作用與凍結管的熱傳導作用達到相對穩(wěn)定，此時溫度場與穩(wěn)態(tài)溫度場非常接近，故可以采用穩(wěn)態(tài)溫度場近似等效瞬態(tài)溫度場。

目前，對于人工凍結穩(wěn)態(tài)溫度場解析解的研究主要集中在無滲流地層。Trupak及Bakholdin分別于1954年及1963年提出單管和單排、雙排凍結管形成的凍結溫度場的解析解；1968年，Sanger通過研究得出單排管凍結溫度場的分布公式；1979年，戸部暢和秋元攻提出多管等間距直線形凍土帷幕穩(wěn)態(tài)溫度場解析解；胡向東等經(jīng)過多年的研究與探索，建立了一套以勢函數(shù)疊加原理為基礎的人工地層凍結穩(wěn)態(tài)溫度場求解方法，并采用該原理，結合一些數(shù)學方法，獲得了一系列溫度場解析解成果，主要有少量管（任意布孔）、排形布管（單排、雙排、三排）、圈形布管（單圈、雙圈）等類別；Wang和王彬等推導了大流速滲透地層單管凍結穩(wěn)態(tài)溫度場解析解，并基于相似模型試驗對解析解的合理性進行了驗證。

目前，對于無滲流地層人工凍結穩(wěn)態(tài)溫度場解析計算的研究成果已較完善，對于滲流場作用下單管凍結穩(wěn)態(tài)溫度場解析解的研究也取得了新的進展，但考慮滲流場影響的多管凍結溫度場解析解尚無相關研究成果。在滲流場作用下，凍結管的吸熱路徑發(fā)生變化，位于水流上游及下游的凍結帷幕的厚度不再相等，該類凍結溫度場解析解與無流速時的對稱凍結帷幕相比存在較大區(qū)別。隨著人工地層凍結法在大流速滲透地層的推廣應用，當前的溫度場解析解已經(jīng)無法滿足實際的工程需要，亟待開展相關研究?；诖?，考慮滲流場影響，對直線單排管非對稱凍結帷幕穩(wěn)態(tài)溫度場的計算公式進行推導；通過水熱耦合物理模型試驗，對公式的合理性進行驗證，并運用得出的公式對不同滲流速度作用下的凍結溫度場的分布規(guī)律進行計算。

1 解析解推導

1.1 基本假設

1）凍結管為直線排列，管壁溫度保持恒定。

2）滲流場的流速、流向、水溫保持不變，水流方向與凍結管連線方向垂直，滲流速度低于凍結帷幕交圈的極限流速。

3）被凍土體為均質材料（凍結范圍內(nèi)土體的熱物理參數(shù)相等）。

4）進入穩(wěn)定凍結階段后，認為凍結溫度場的變化速率較慢，在某一時刻的凍結溫度場保持不變。

5）凍結帷幕向上游、下游及兩側方向的擴展半徑分別為，均為滲流速度

v

的函數(shù)；對應同一種流速

v

，進入穩(wěn)定凍結階段后的及保持不變，故在公式推導過程中可忽略流速

v

。

6）直線單排凍結管最外側的凍結鋒面等效為兩個軸長不同的1/4橢圓，如圖1所示。

圖1 滲流場作用下3管凍結穩(wěn)態(tài)溫度場計算模型Fig. 1 Calculation model of steady-state freezing temperature field of three pipes under the action of seepage field

該部分凍結鋒面的軌跡方程表示為：

7）由于滲流場的作用，下游凍結帷幕的厚度不均勻，其中，最大厚度為ξ，最小厚度為ξ。為便于計算，假設下游凍結帷幕的厚度相等ξ=(ξ+ξ)/2；上游凍結壁的厚度為ξ。

1.2 公式推導

不同流速條件下的溫度場解析解形式相同，滲流速度通過影響邊界條件（凍結帷幕的擴展范圍）來影響公式的計算結果，因此在公式推導過程中可忽略流速。

1）凍結帷幕中間部分溫度場解析解

將滲流場作用下直線單排凍結管形成的凍結帷幕的中間部分等效為直線單排非對稱發(fā)展的凍土帷幕，溫度場數(shù)學模型為：

式中，

T

為凍結溫度，

T

為凍結管壁溫度，

r

為凍結管半徑，

i

為凍結管的數(shù)量。

式（4）～（5）中，

T

及

T

為待定常數(shù)。

通過求解得出：

2）凍結帷幕兩側部分溫度場解析解

單管凍結穩(wěn)態(tài)溫度場的方程為：

當?shù)貙又袩o滲流場作用時，單管凍結鋒面的形狀為圓形，

L

為定值；當?shù)貙又写嬖跐B流場時，單管凍結鋒面的形狀不再是圓形，

L

隨著流速及位置而改變，故需要對

L

的表達式進行推導。經(jīng)過凍結區(qū)域內(nèi)一點

M

(

x

,

y

)及最外側凍結管（0,

d

）的直線方程為：

聯(lián)立式（1）、（10）可以得出經(jīng)過

M

(

x

,

y

)的直線與凍結鋒面的交點為：

位于下游區(qū)域（

x

＞0）凍結鋒面上的點到凍結管的距離可以表示為：

凍結區(qū)域內(nèi)的點

M

(

x

,

y

)到凍結管的距離表示為：

同理，可得到上游區(qū)域（

x

＜0）的凍結鋒面到凍結管的距離為：

凍結區(qū)域內(nèi)點

N

(

x

,

y

)到凍結管的距離表示為：

由于凍結鋒面上的任一點，同樣可以視為半徑為

L

的圓上的一點，將

L

及

L

分別代入式（9），考慮點

M

(

x

,

y

)及

N

(

x

,

y

)選取的任意性，得出對應凍結區(qū)域的穩(wěn)態(tài)凍結溫度場的解析表達式為：

當ξ=ξ=ξ，且

d

=0時，式（16）退化為單管凍結溫度場的Trupak公式。

在式（8）的推導過程中，假定凍結壁的厚度與凍結管間距之間的比值無限大，但在實際工程中，尤其是試驗中，無法滿足該條件；而在式（16）的推導中，沒有對凍結帷幕的厚度進行限定，因此可以基于式（16）的計算結果對式（8）進行修正。

當

y

=

d

時：

式（8）的計算結果為：

式（16）的計算結果為：

對式（8）進行修正，兩種修正系數(shù)的形式為：

對兩種修正系數(shù)的修正效果進行計算分析后發(fā)現(xiàn)：采用修正系數(shù) γ可以保證凍結鋒面處的溫度計算值與凍結溫度

T

相等，但當下游凍結壁厚度遠大于上游凍結壁厚度時，若采用 γ進行修正，位于兩根凍結管中間位置的溫度計算數(shù)據(jù)會出現(xiàn)突變點，溫度計算結果誤差較大；采用 γ進行修正時，溫度變化曲線較為平滑，但是不能保證凍結鋒面處的溫度與凍結溫度

T

相等，并且上游最低溫度計算值并不出現(xiàn)在凍結管管壁處。因此，在上游區(qū)域（

x

＜0），采用修正系數(shù) γ進行修正；在下游區(qū)域（

x

＞0），采用修正系數(shù)γ進行修正。

修正后，滲流場作用下直線單排管非對稱凍結帷幕穩(wěn)態(tài)溫度場中間部分的溫度計算公式為：

3）凍結帷幕厚度計算公式

實際工程中，經(jīng)常需要根據(jù)測點的溫度判斷凍結帷幕的厚度，若在凍結壁的上、下游各布置一個測溫點，測溫點坐標分別為（

X

,

Y

）、（

X

,

Y

），對應的測點溫度分別為

T

及

T

，對式（21）進行變換，可以得出凍結壁厚度的計算公式為：

4）凍結帷幕平均溫度計算公式

平均溫度是凍結帷幕的重要參數(shù)。若凍結孔的間距相等，凍結帷幕的平均溫度等于凍結帷幕區(qū)間0≤

y

≤0.5

d

上的平均溫度。因此，凍結帷幕的平均溫度計算公式可以表示為：

式（23）的積分過程較為復雜，在實際工程中可采用等效梯形法計算直線單排管凍結帷幕的平均溫度。由于在實際工程中凍結管的數(shù)量較多，因此，在平均溫度的計算過程中，可以忽略凍結帷幕兩側形狀不規(guī)則對計算結果的影響。目前，已有研究成果證明距離凍結管

d

/5處的截面平均溫度可近似等效為整個凍結帷幕的平均溫度，則計算公式可簡化為：

式中，

T

表示凍結帷幕截面與軸面交點處的溫度。

2 解析解合理性驗證

2.1 模型試驗設計

為對解析解的合理性進行驗證，基于推導過程中的基本假設，設計了物理模型試驗系統(tǒng)，如圖2所示。該試驗系統(tǒng)由多孔介質試驗區(qū)、凍結系統(tǒng)、滲流場模擬系統(tǒng)、數(shù)據(jù)采集系統(tǒng)等組成。

圖2 試驗系統(tǒng)3D示意圖Fig. 2 3D schematic diagram of test system

為滿足基本假設1），設計箱體的尺寸為2 500 mm×2 000 mm×1 000 mm，兩端各安裝8個規(guī)格為40 mm×2 mm的鋼管，分別作為進水管及出水管。在箱體中軸線處按照直線布置3根外直徑為42 mm的凍結管，管間距為400 mm。試驗中，保持制冷機組的酒精溫度為–32℃，流量為2.5 m/h。

為滿足基本假設2），兩道100目的濾網(wǎng)將設計箱體沿長度方向分成3部分，中間部分為主體實驗室，其長度為2 000 mm，用來容納多孔介質。在箱體的兩端各設置一段250 mm的緩沖室，緩沖室內(nèi)填滿石子，從水管中流出的水流會通過緩沖層的縫隙流向中間區(qū)域的多孔介質，以保證水流能夠沿著箱體截面均勻進入中間的主體試驗區(qū)。通過恒壓變頻泵控制水的流量，保證每次試驗過程中的水流量保持不變。儲水箱為恒溫水箱，從而保證整個循環(huán)過程中水流的溫度不變。

研究共進行4組試驗，試驗中控制的滲流速度分別為0、3、6、9 m/d，對應的流量及水溫見表1。

表1 模型試驗中水流參數(shù)的設置值
Tab. 1 Setting values of water flow parameters in the model test

編號滲流速度/(m·d–1)流量/(m3·h–1)水溫/℃1 0 0—0.2515 3 6 0.5015 2 3 4 9 0.7515

為滿足基本假設3），選用粒徑為（1±0.15）mm的均質圓粒砂作為多孔介質模擬材料（圖3），砂的粒徑均勻、無雜質；在填充過程中，每填充5 cm壓實一次，可近似認為多孔介質為均質材料。

圖3 多孔介質模擬材料Fig. 3 Similar material of porous medium

試驗系統(tǒng)的主測試平面位于箱體500 mm層位，分為A、B、C、D、E、F、G共7條軸線，每條軸線分布有13個測點，測點間距為50 mm，共計91個測點；輔助測試平面位于箱體400 mm層位，分為H、I、J共3條軸線，每條軸線上有27個測點，測點間距為50 mm，共計81個測點。測點布置情況如圖4所示。

圖4 溫度測點布置示意圖Fig. 4 Schematic diagram of temperature measurement points arrangement

多孔介質與箱體邊界的接觸面會形成邊界效應，并隨著試驗次數(shù)的增加，多孔介質可能會與箱體上蓋之間形成縫隙，水流會向阻力小的方向流動，最終導致多孔介質中的滲流速度與試驗設計值存在較大偏差。為消除該效應對滲流場帶來的影響，在箱體的底面、側面及上蓋上設置直徑為25 mm的阻水條（底板及上蓋各設置2根，兩個側面各設置3根），防止水流沿著箱體邊界流動，如圖5所示。

圖5 箱體邊界阻水條設置情況Fig. 5 Setting of water bar at box boundary

為防止外界的環(huán)境溫度對凍結過程產(chǎn)生干擾，對箱體外表面及凍結管路進行保溫處理。首先，緊貼箱體外表面布置一層30 mm的橡塑保溫板；隨后，采用40 mm的聚氨酯保溫板將整個箱體包裹在其中，接縫處用氯丁膠緊密粘合，外露的凍結管及凍結干管表面包裹一層30 mm的橡塑保溫層，保溫層的布置情況如圖6所示。

圖6 箱體保溫層設置情況Fig. 6 Setting of insulation layer of box

整體模型試驗系統(tǒng)如圖7所示。

圖7 水熱耦合模型試驗系統(tǒng)Fig. 7 Hydrothermal coupling model test system

2.2 模型試驗數(shù)據(jù)分析

通過模型試驗得出凍結帷幕的交圈時間如表2所示（以兩管中間測點溫度到達–1℃為交圈判斷依據(jù)）。分析后發(fā)現(xiàn)，凍結帷幕的交圈時間與滲流速度之間存在近似指數(shù)函數(shù)關系。采用指數(shù)函數(shù)對試驗結果進行擬合，結果如圖8所示。

表2 不同流速條件下凍結帷幕的交圈時間
Tab. 2 Closure time of frozen wall under different seepage velocity

滲流速度/(m·d-1)交圈時間/min延遲時間/min 0 720—3 80585 6 1 760955 9 4 3002 540

圖8 凍結帷幕交圈時間的擬合曲線Fig. 8 Fitting curve of closure time of frozen wall

基于擬合結果，提出滲流場作用下凍結帷幕交圈時間的預測公式為：

式中：

v

為滲流速度；

W

、

t

、

m

為擬合參數(shù)，在本文試驗中，分別為120.46、526.57、–2.62。定義凍結帷幕向下游的最大擴展半徑

R

與凍結帷幕向上游的最大擴展半徑

R

的比值

R

/

R

為凍結帷幕的非對稱系數(shù)，其隨著滲流速度的變化規(guī)律如圖9所示。

圖9 凍結帷幕非對稱系數(shù)隨著滲流速度的變化規(guī)律Fig. 9 Variation of asymmetry coefficient of frozen wall with seepage velocity

由圖9可知，隨著滲流速度的增加，凍結帷幕的非對稱性逐漸加劇。

2.3 公式計算結果與試驗結果對比

為驗證本文推導的解析解的準確性，將進入穩(wěn)定階段后的B、C軸及I軸（圖4）的溫度計算結果與模型試驗結果進行對比，如圖10～12所示。

圖10 B軸線上測點溫度的計算結果與試驗結果對比Fig. 10 Comparison of calculation results and test results of temperature at the measurement points on the B axis

通過對比發(fā)現(xiàn)，B軸上測點溫度的計算結果與試驗結果吻合度較高，僅少數(shù)測點的誤差大于1℃。I軸線上位于凍結管外側（等效為兩個1/4橢圓部分）的測點的溫度計算結果與試驗結果的吻合度較高，但該軸線上位于凍結管中間部分測點的計算結果誤差較大，個別測點的誤差超過2℃。C軸線上測點溫度計算結果的整體分布規(guī)律與試驗結果一致，但溫度數(shù)值存在一定誤差，其中，下游區(qū)域計算結果的誤差大于上游區(qū)域，且隨著水流速度的增加，計算結果的誤差增大。

公式計算結果出現(xiàn)誤差的原因如下：位于整個B軸和I軸外側的測點主要受單根凍結管的熱傳導和地層中水流對流傳熱的影響，而整個C軸和I軸中間部分測點受多根凍結管熱傳導疊加和水流對流傳熱的影響；相鄰凍結管中間區(qū)域砂層凍結過程更加復雜，穩(wěn)態(tài)溫度場計算公式無法同時考慮多管熱傳導疊加及水流對流傳熱作用，因此兩管中間區(qū)域的計算結果存在一定的誤差。

3 滲流場作用下3管凍結穩(wěn)態(tài)溫度場計算

已知凍結帷幕向上、下游的擴展厚度及凍結管管壁的溫度，可求得不同流速條件下直線單排3管穩(wěn)態(tài)凍結溫度場的分布規(guī)律（圖13），以及凍結帷幕各截面的溫度分布規(guī)律（圖14）。

圖13 不同流速條件下3管凍結穩(wěn)態(tài)溫度場計算結果Fig. 13 Calculation results of steady-state temperature field of three-pipe freezing under different flow rate conditions

圖14 沿著水流方向不同截面溫度分布規(guī)律Fig. 14 Temperature distribution law of different cross-sections along the flow direction

圖11 I軸線上測點溫度的計算結果與試驗結果對比Fig. 11 Comparison of calculation results and test results of temperature at the measurement points on the I axis

圖12 C軸線上測點溫度的計算結果與試驗結果對比Fig. 12 Comparison of calculation results and test results of temperature at the measurement points on the C axis

由圖13及14可知：凍結帷幕各截面的溫度隨著與凍結管距離的增加而升高，越靠近凍結壁中間位置，各個截面的溫差越大，除凍結管所在截面（

y

=0）之外，其他截面的溫度分布曲線都類似于梯形。隨著水流速度的增加，凍結壁的不對稱性加劇，但截面上溫度變化規(guī)律較為接近，距離凍結管所在軸線的距離越遠，溫度曲線趨于平緩，即距離凍結管越遠，溫降速率越低。

根據(jù)凍結帷幕平均溫度計算公式，得出試驗中不同流速條件下凍結帷幕的平均溫度，見表3。

表3 凍結帷幕平均溫度計算結果
Tab. 3 Calculation results of average temperature of frozen wall

流速/(m·d–1)凍結時間/h厚度/mm平均溫度/℃0 32540–10.10 3 44607–10.54 6 48559–10.22 9 96551–10.00

凍結試驗中，由于流速越大，對應凍結帷幕的交圈時間越長，因此，每個流速對應的凍結試驗持續(xù)時間并不相同，其中：無流速時試驗持續(xù)時間為32 h；流速為9 m/d時，試驗持續(xù)時間為96 h。當?shù)貙又袩o滲流場時，凍結32 h后，凍結帷幕的平均溫度為–10.10 ℃。當?shù)貙又写嬖跐B流場時，由于水流的對流傳熱作用，凍結管的傳熱路徑發(fā)生了變化，大量冷量在下游集聚，此時凍結帷幕的厚度受到滲流速度及凍結時間的影響，凍結帷幕的平均溫度變化過程較為復雜，但整體仍然表現(xiàn)出隨著凍結壁厚度的增加而降低。當滲流速度為3、6、9 m/d時，凍結試驗持續(xù)時間分別為44、48、96 h，對應的凍結壁厚度分別為607、559、551 mm，通過計算得出凍結壁的平均溫度分別為–10.54、–10.22、–10.00 ℃。

4 結 論

1）基于穩(wěn)態(tài)溫度場的求解理論，采用分段等效的方法，推導得出定向滲流誘導的非對稱凍結帷幕穩(wěn)態(tài)溫度場解析解，并提出凍結帷幕厚度及平均溫度的計算方法。

2）采用自主構建的水熱耦合物理模型試驗系統(tǒng)，對不同流速條件下直線單排管凍結溫度場的形成規(guī)律進行了研究。結果表明，凍結帷幕的交圈時間及非對稱系數(shù)隨著流速的增加急劇增大。

3）基于模型試驗的結果對解析解的合理性進行驗證，對比發(fā)現(xiàn)，公式計算結果與模型試驗結果整體的吻合度較高，但由于多根凍結管熱傳導及水流對流傳熱的疊加作用，位于兩管中間區(qū)域的溫度場計算結果存在一定誤差。

4）采用本文推導的公式對滲流場作用下凍結溫度場的分布規(guī)律進行了計算分析，結果表明：當?shù)貙又写嬖跐B流場時，凍結溫度場變化過程較為復雜，但平均溫度整體上仍表現(xiàn)出隨著凍結帷幕厚度增加而降低的規(guī)律。