丁法龍，茅澤育

(清華大學水利水電工程系，北京 100084)

物體或結構受熱之后，由于受到外部或內部的約束而不能實現自由膨脹，物體或結構內部將產生溫度應力，又稱熱應力。在寒冷地區(qū)，水庫、河流、湖泊、水塘等水體在冬季都有冰凍現象發(fā)生并形成一定厚度的冰層(淡水冰)，冰具有固態(tài)物體熱脹冷縮的物理性質，冰在發(fā)生膨脹時，會受到周圍或內部結構物(如水庫邊坡護岸、壩體、閘門、海域港口防波堤等)不同程度的約束從而在冰層內部產生應力，冰層反過來對結構物施以擠壓破壞作用，即冰壓力[1]。需要說明的是，由水凍結成冰的相變過程中，盡管體積也會發(fā)生膨脹，但此時冰的材料強度處于較低的水平上，不會對結構產生明顯的影響[2]；然而已有大量觀測資料表明，冰在升溫過程中會產生較大的溫度應力，該溫度應力引起的冰壓力可導致結構物不同程度的破壞，是寒區(qū)水工建筑物的關鍵設計參數[3]。本文主要研究由冰溫升高引起的冰層溫度應力。

冰的物理特性十分復雜，研究表明冰溫度應力受多種因素的影響，包括冰層的熱力狀況、冰層的約束條件、冰層厚度等，這使得冰溫度應力的估算存在諸多不確定因素。我國海冰物理力學性質[4]和動冰力[5?7]的研究相對較為系統，但針對淡水冰尤其是靜冰力的研究仍然有很大缺口，國內對該問題的研究基本上采用的是現場量測冰壓力的大小，經過數據分析提出相應的經驗公式，將冰壓力表達為冰厚和結構物幾何參數的函數[8]，這些經驗公式具有一定的工程應用基礎，但它們的適用范圍小，誤差大，而且公式中一些經驗參數在實際中難以確定，是一種粗略的冰力估算方法。通過建立冰層應力-應變的力學模型(本構模型)，對冰溫度膨脹力的作用機理進行分析，是一類更為深入的研究思路。史慶增等[9]將層合板理論引入冰層溫度應力的分析求解中，提出了一種求解對邊約束冰層溫度膨脹力的計算方法，并對海冰進行了計算；黃焱等[10]采用有限元方法分析了單層淡水冰的溫度膨脹力，并對水工建筑物中常見的形式進行具體計算。這些基于應力-應變力學模型的冰力數值計算研究中均假定了冰的彈性本構關系。實際上，冰的本構關系與應變率明顯相關，靜冰溫度應力發(fā)生時，彈性變形階段只有很短的時間，因此，利用彈性本構關系求解靜冰溫度應力問題可能導致誤差。黃文峰等[11]以紅旗泡水庫為例，開展不同方向上冰表面形變的原位觀測，嘗試改進現有的冰溫度應力計算模型來預測表層應力。國外靠近北極地區(qū)的國家(瑞典、挪威、加拿大等)，對湖泊及大型水庫的冰溫度應力計算開展了大量的試驗和理論研究，并采用了多種流變學本構模型[12?16]用于冰溫度應力計算，其中大多是將兩個基本的粘彈性本構模型—Maxwell模型和Kelvin模型以串聯或并聯的方式進行組合，然后，根據冰力原型觀測數據確定模型中的待定參數，用于當地的靜冰溫度應力計算，是至今較常采用的靜冰力計算方法。但是這些冰力模型中的參數取值在很大程度上依賴于數據獲取的特定環(huán)境條件，不能直接應用于中國北方寒冷地區(qū)的靜冰溫度應力計算。

此外，迄今為止，國內外冰溫度應力的研究多集中于湖冰和庫冰，關于淺水湖塘中冰層溫度應力的研究極少見于報道。與大型深水湖庫相比，小型淺水湖塘的傳熱條件不同，從而導致冰層熱力狀況有所不同；另外，小型湖塘的面積更小，周邊圍岸對冰層的約束作用更為強烈。因此，小型淺水湖塘中冰層溫度應力呈現特有的時空變化規(guī)律。

綜上，國外的靜冰力模型不可直接應用于中國東北地區(qū)的靜冰力計算，有關塘冰溫度應力的研究更是極為少見。針對上述2個研究背景，為探究我國東北地區(qū)小型靜態(tài)水域的冰層溫度應力變化規(guī)律，本文采用原型觀測研究手段，結合理論分析及計算，以黑龍江省大慶市青花湖6號水塘為研究對象，對水塘冰層內部溫度和應力進行了觀測及分析，基于冰的粘彈性本構關系提出了一種塘冰溫度應力計算模型，并結合實測數據回歸得到了模型中的待定參數，該模型計算值與實測值吻合良好，可為寒區(qū)小型淺水湖塘及小型水庫的冰溫度應力計算提供參考。

1 溫度應力原型觀測

1.1 觀測場地概況

冰情原型觀測現場位于黑龍江省大慶市肇源縣青花湖。青花湖，原名庫里泡，位于大慶市大同區(qū)、肇州縣、肇源縣三地交界處，地處124.78°E～124.87°E，45.77°N～45.91°N，是安肇新河上的最終控制水體，湖區(qū)從南到北由幾個相連的湖盆組成，南北長約16 km，東西最寬處約5 km，湖區(qū)水面面積64.05 km2，平均水深1.2 m～1.5 m，最大水深約2.8 m。青花湖水庫每年10月下旬至11月初開始結冰，冰層到次年4月中下旬融化，冰期5～6個月，最大結冰厚度1.0 m～1.2 m。每年7月～8月份給湖內補水，冬季結冰期停止補水，基本沒有水位變動，因此青花湖的冬季自然環(huán)境接近靜態(tài)水體，提供了靜冰觀測的理想條件。

青花湖西端的部分水域被人為地以漿砌石擋墻劃分為若干水塘進行漁業(yè)養(yǎng)殖，本次觀測選取形狀規(guī)整、塘底高程起伏不大的青花湖6號水塘為研究目標。冰溫度應力主要發(fā)生在連日持續(xù)升溫的條件下，根據對當地以往多年氣溫和冰情資料的分析，選擇觀測時段為2018年2月12日?3月14日，對當地氣溫、冰溫及冰溫度應力進行了持續(xù)觀測。

青花湖水庫6號水塘呈東西長度120 m、南北寬度82 m的長方形布置，平面布置如圖1(a)所示；四周布有高出水面0.8 m左右的漿砌石擋墻，可視為對冰層的垂直剛性約束邊界。

圖1 原型觀測布置圖Fig.1 Sketch of prototype observation layout

1.2 測量裝置與方法

冰溫及冰溫度應力的測量裝置包括：1)冰溫測量采用羅漢姆公司生產的RH-8068熱電阻pt100鉑熱電阻溫度傳感器，測溫精度±0.1 ℃；2)冰壓力測量采用Interlink Electronics 公司生產的FSR402電阻式壓力傳感器，這款壓力傳感器是將施加在FSR傳感器薄膜區(qū)域的壓力轉換成電阻值的變化，從而獲得壓強信息，精度為1 kPa；3)氣溫數據由通威集團青花湖光伏站的A753 WS無線自動氣象站提供；4)冰溫及冰壓力測量數據采集至CR1000X數據采集儀，數據采集時間間隔均為10 min/次。為探究冰溫度應力隨邊界約束條件的變化，在水塘內設置了4個不同的平面測點位置，即自水塘東西擋墻的中垂線，距離北側擋墻2 m，自西向東每隔15 m設置一個觀測位置，分別記為測點1、測點2、測點3和測點4(如圖1(a)所示)。

在每個觀測點處，從冰面至冰下分別布設一條冰溫和冰溫度應力觀測鏈(如圖1(b)所示)，觀測鏈上安裝有16組pt100溫度傳感器和FSR402電阻式壓力傳感器，觀測鏈采用木質框架固定，選擇木質材料的原因是木頭的導熱系數低，對冰的熱影響較小。溫度傳感器探頭和壓力傳感器探頭總是成對地處于同一冰深處，以自動記錄該垂向位置的冰溫和應力變化?？紤]到冰層上部的溫度變化速率和垂向梯度大于下部，因此冰層上部的溫度傳感器和壓力傳感器布置較為密集，每5 cm布置一組，下部每隔10 cm布置一組。按照最大冰厚度1.2 m設計，傳感器組在冰面以下的具體布設位置為：5 cm、10 cm、15 cm、20 cm、25 cm、30 cm、35 cm、40 cm、50 cm、60 cm、70 cm、80 cm、90 cm、100 cm、110 cm、120 cm。

觀測鏈安裝前，先在冰層上沿平行于擋墻方向鋸開一條貫穿的窄縫，將已經固定好溫度傳感器探頭和壓力傳感器探頭的觀測鏈豎直放入窄縫之中，規(guī)定浸水后每個壓力傳感器的讀數為該壓力傳感器的零點，觀測鏈導線連接在木質支架上，木質支架固定于冰層上。冰溫和冰壓力觀測鏈于2018年2月2日在青花湖水庫6號水塘安裝完成，當日水塘冰厚為97 cm。觀測鏈的溫度探頭和壓力探頭隨著窄縫內部冰生長被凍結到冰內，2018年2月12日窄縫凍結密實后開始進行數據采集，以保證傳感器探頭與冰充分接觸，從而精確感知該深度處冰溫及冰應力的變化。冰厚觀測采用冰鉆人工穿孔后用塔尺直接量測。

2 理論模型

2.1 冰的熱變形

預測冰力大小的精確性很大程度上取決于冰力學模型的合理性，從冰的溫度應力產生機理分析，當冰升溫產生的膨脹變形受到周界約束時，冰層內部產生溫度應力，溫度應力的大小應恰好抵消這部分原本應該產生的溫度應變，因此，冰溫度應力計算主要取決于冰的熱變形關系和靜冰應力-應變本構關系。

由于熱膨脹只產生線應變，為剪切應變?yōu)榱?，對于三維問題，冰層溫度升高引起冰的熱變形關系可表示為：

式中： ε為溫度應變； ΔT/ (℃)為溫度增量；α為冰的溫度膨脹系數，可取α=(54+0.18T)×10?6℃?1。

如果認為膨脹系數為常數，對式(1)兩邊對時間求導得：

江西茶葉是以綠茶為主導，綠茶的總產量占據總產茶量的80%以上，而紅茶的總產量僅占據10%左右，而其他種類的茶葉占比更是不足10%。這就在一定程度上限制了江西以綠茶出口為主的發(fā)展空間。

由式(2)可見，冰溫變化率與應變速率為線性關系，因此，冰溫變化引起的溫度應力求解又歸結為確定對應的應變速率下的冰本構關系。

2.2 流變性本構關系

靜冰溫度應力發(fā)生在低應變速率下(10?7s?1～10?8s?1)，此 時 冰 表 現 為 典 型 的 粘 彈 性(即 流 變性)，在眾多的靜冰流變學本構模型中，瑞典學者Bergdahl[12]基于Maxwell單元提出了一種改進模型，該模型僅有3個待定參數且能較好地描述冰的力學特性，是最簡單的模型之一，其可靠性已得到多項研究[13,15?16]的實測驗證。在Bergdalh模型中，冰應力-應變關系被描述成由表征彈性的彈簧元件和表征粘性的非線性阻尼器元件組成的串聯模型，其表達式為：

聯立式(2)、式(3)可得：

式中： σ/kPa為某點的溫度應力；t/s為時間；T/(℃)為冰溫；D/(m2/s)為蠕變速率；K/m?2為粘性變形系數；n為應力指數； σ0是為了實現應力無量綱化引入的常應力值，取為100 kPa。Bergdahl[12]和Cox[13]等基于?10 ℃條件下的S1柱狀冰來確定冰的力學參數K、D、n，本研究中采用了他們的取值。E是冰的彈性模量，可描述為E=E0(1?0.012T)，其中Bergdahl[12]取E0=6.1 GPa，Cox[13]則取E0=4 GPa。

Bergdahl[12]在模型推導過程中借用了阿倫尼烏斯公式，即化學反應速率常數隨溫度變化關系的經驗公式，來比擬冰的蠕變速率隨溫度的變化關系。Cox[13]認為化學反應與冰的蠕變過程區(qū)別較大，這種比擬缺乏一定的解釋力，并通過實測數據回歸，以相對溫度的指數函數形式替換Bergdahl模型中的KD，即令：

式中：T1=?1 ℃；m=1.92[13]；β均為擬合參數。

由以上推導過程可以證明，參數A即是冰彈性模量E與冰熱膨脹系數 α的乘積，即：

由于A是由兩個具有明確物理意義的參量組合而成，故在對式(6)擬合分析時，不再把參數A作為純經驗參數進行擬合確定，而是以某個常量代替，以避免回歸分析中A、B兩個參數的的相互影響。E和 α都是隨冰溫T變化的量，根據溫度場觀測結果，冰溫T集中分布于?15 ℃～0 ℃區(qū)間內，取E0=4 GPa，由式(7)計算得到該溫度區(qū)間內變量A的平均值(積分中值)為229 kPa/℃，本研究中以該常量代替參數A，主要分析模型中待定參數B的取值。

另外需要說明的是，以上分析過程均是基于冰內應力為正向的假設，由于壓力傳感器的零點設置為浸水時的讀數，冰在升溫膨脹時內部產生壓應力，此時壓力傳感器測得正值，但降溫收縮時壓力傳感器有可能測得負值，為便于數據采集和分析，在式(6)中引入符號函數，得到本文所采用的冰溫度應力計算模型：

式(6)、式(8)中：A=229 kPa/℃；B為待定參數，單 位 為kPa/d；m=1.92[13]；n=3.651[13]；σ0=100kPa ；T1=?1 ℃； sign(σ) 為 σ的符號函數，當σ>0 時取1， σ<0時取?1；其他符號意義同上文式(1)～式(7)。

冰溫和應力觀測的數據采集時間間隔為10 min，式(8)可采用顯式時間積分算法由1stOpt商業(yè)數學軟件求解，參數B由單純形算法回歸得到。

3 結果與分析

氣溫是影響冰層生消、冰層溫度場及應力場的最主要熱力因素，冰溫及應力分布與氣溫變化緊密相關。氣溫觀測資料(圖2)表明，該觀測區(qū)域于2018年2月初開始出現日平均氣溫回升，本研究進行數據采集的時間段為2018年2月12日?3月14日為比較穩(wěn)定的連續(xù)升溫階段，符合產生冰溫度應力的氣溫條件，觀測期內冰層厚度穩(wěn)定在97 cm～102 cm。下面根據實際觀測數據，分別分析冰層溫度和溫度應力的時空分布規(guī)律。

圖 2 2017年11月?2018年4月觀測地日平均氣溫變化Fig.2 Daily mean air temperature from Nov. 2017 to Apr. 2018

3.1 冰層溫度

冰層溫度是冰的基本狀態(tài)參量，導致溫度應力產生的前提條件是溫度變化與約束作用，冰溫在冰層內部的變化直接影響到冰層溫度應力的大小和分布。圖3為2018年3月1日觀測點1處不同深度冰層的冰溫及氣溫變化，由圖可見，冰溫變化主要取決于上部氣溫的波動情況，各層冰溫隨時間的波動較氣溫更為平緩，表層(5 cm)冰溫變化隨氣溫變化的響應最為敏感，且隨著冰層深度的增加，冰溫波動幅度逐漸減小，冰溫也逐步升高。在冰層深度30 cm以內，冰溫呈現較為明顯的日周期變化，在30 cm以下，冰溫的日變化特征逐漸減弱，相對于氣溫的滯后時間更長。冰溫與氣溫最小值出現在6時～7時，冰溫比氣溫大約滯后1 h，最大值出現在下午12時～17時，冰溫滯后約2 h。

圖3 2018年3月1日觀測點1處不同深度冰層冰溫與氣溫日變化Fig.3 Diurnal variation of ice temperature of different depths and air temperature at station 1 on March 1, 2018

冰層表層溫度由大氣和冰面之間的熱交換所決定，并與冰內的熱傳導過程相耦合，同時還受到太陽輻射、云遮、風速、覆雪等因素的的影響，但有研究[7]表明，冰面冰溫與氣溫的相關性顯著，通過統計觀測期內水塘冰層表層(冰面以下5 cm)冰溫及同時刻氣溫數據，繪成圖4，由圖可見，每日不同時刻水塘冰層表層冰溫均與氣溫呈現良好的線性關系，且變化率在0.38～0.56。

圖4 冰表面以下5 cm處的冰溫與氣溫關系Fig.4 The relationship between ice temperature at 5 cm depth and air temperature

圖5 觀測點3處2018年2月21日不同時刻的冰溫垂向分布Fig.5 Vertical distribution of ice temperature at different times at station 3 on Feb 21, 2018

3.2 冰溫度應力

采用文獻[13]中的的取值n=3.651，m=1.92及A=229 kPa/℃，以應力觀測值對式(8)中的待定參數B進行擬合回歸，得到不同測點和不同冰深位置處的冰溫度應力經驗計算模型。

圖6為觀測期內不同測點位置處冰表層(冰深5 cm處)溫度應力觀測值與模型計算值，由圖可見最大的應力觀測值出現在2018年2月19日的測點4處，為322 kPa。通過對比應力觀測值與計算值發(fā)現，盡管式(8)形式簡單且只用到一個擬合參數B，但模型計算結果與觀測結果吻合良好，在氣溫穩(wěn)定回升的2018年2月20日?2018年3月10日期間，測點1、測點2、測點3、測點4處的冰溫度應力模型計算值與實測值之間的平均相對誤差分別為?4.63%、?5.19%、?7.26%、?5.37%，均在±10%以內，表明該模型用于計算冰層溫度應力具有較好的精度。

圖6 不同測點位置在冰深5 cm處的靜冰溫度應力觀測值與計算值對照Fig.6 Observed and calculated ice stress at 5 cm deep at different stations

與Cox[13]的研究結果(B=27 kPa/d)相比，本文的擬合參數B值較大且明顯依賴于測點位置。即從水塘側面(測點4，B=8646 kPa/d)向水塘中心(測點1，B=469 kPa/d)，隨著邊界約束作用的減弱，B值遞減。根據2.2節(jié)理論模型部分的介紹，可以預期參數B僅是冰的固有材料性質，然而，觀測數據表明B與測點位置相關，限于本試驗測點位置布置的單一性，尚無法通過觀測數據確定B的成因，關于B的具體影響因素及分布規(guī)律可通過更加全面深入的觀測數據進行分析。

圖7為2018年2月21日不同時刻，冰深10 cm處的冰應力在不同測點處的分布，可見從測點1到測點4，整體應力水平呈遞增趨勢，這是因為兩邊擋墻對水塘冰層的邊界約束作用越來越強所致。

圖7 2018年2月21日冰深10 cm處冰溫度應力水平分布Fig.7 Horizontal distribution of ice stresses at 10 cm depth at different times on February 21, 2018

圖8為測點3處，2018年2月21日不同時刻的冰應力垂向分布廓線，由于各層冰溫變化值及冰層狀態(tài)不同，冰層不同深度處的溫度應力也不同。冰層表面為自由面，受到的約束較小，因而產生的應力稍小，隨著深度的增加冰層溫度應力增大，并在冰深10 cm～30 cm處達到最大，其下又隨深度逐漸減小，綜合觀測期內的所有冰應力垂向分布廓線，冰應力只在0.7 m以內的冰層上部產生，這是因為冰層內部的熱傳導主要受氣溫主導，冰層上部與大氣的熱交換較下部強烈，冰溫變化率沿冰深逐漸衰減，由式(2)可知當冰溫變化率為零時，冰的應變速率為零，從而不會產生溫度應力。

圖8 測點3處，2018年2月21日不同時刻冰溫度應力垂向分布Fig.8 Vertical distribution of ice stress at station 3 and different times on February 21, 2018

4 結論

采用原型觀測研究手段，以黑龍江省大慶市肇源縣青花湖6號水塘為研究場地，對水塘冰層內部溫度場和應力場進行了研究，主要結果及結論如下：

(1)冰溫變化主要取決于上部氣溫的波動情況，表層(5 cm)冰溫變化隨氣溫變化的響應最為敏感，且隨著冰層深度的增加，冰溫波動幅度逐漸減小，冰溫也逐步升高。在冰層深度30 cm以內，冰溫呈現較為明顯的日周期變化，在30 cm以下，冰溫的日變化特征逐漸減弱。

(2)表層冰溫均與氣溫呈現良好的線性關系，且變化率在0.38～0.56；冰溫垂向分布廓線的上半部分(30 cm以上)冰溫隨時間呈現復雜的變化，下半部分(30 cm以下)冰溫沿冰深基本呈線性分布。

(3)基于Bergdahl冰流變性本構關系提出了一種冰溫度應力計算模型，并結合實測數據給出了模型中參數的確定方法，擬合參數B值明顯依賴于測點位置；盡管模型形式簡單且只用到一個擬合參數，但模型計算結果與觀測結果吻合良好。

(4)冰層上部產生的溫度應力比下部大，冰層表面為自由面，受到的約束較小，因而產生的應力稍小，隨著深度的增加冰層溫度應力增大，并在冰深10 cm～30 cm處達到最大，其下又隨深度逐漸減小，綜合觀測期內的所有冰應力垂向分布，冰應力只在0.7 m以內的冰層上部產生；從測點1～測點4，整體應力水平呈遞增趨勢。