寒区水塘冰温度应力原型观测及分析

丁法龙，茅泽育

(清华大学水利水电工程系，北京 100084)

物体或结构受热之后，由于受到外部或内部的约束而不能实现自由膨胀，物体或结构内部将产生温度应力，又称热应力。在寒冷地区，水库、河流、湖泊、水塘等水体在冬季都有冰冻现象发生并形成一定厚度的冰层(淡水冰)，冰具有固态物体热胀冷缩的物理性质，冰在发生膨胀时，会受到周围或内部结构物(如水库边坡护岸、坝体、闸门、海域港口防波堤等)不同程度的约束从而在冰层内部产生应力，冰层反过来对结构物施以挤压破坏作用，即冰压力[1]。需要说明的是，由水冻结成冰的相变过程中，尽管体积也会发生膨胀，但此时冰的材料强度处于较低的水平上，不会对结构产生明显的影响[2]；然而已有大量观测资料表明，冰在升温过程中会产生较大的温度应力，该温度应力引起的冰压力可导致结构物不同程度的破坏，是寒区水工建筑物的关键设计参数[3]。本文主要研究由冰温升高引起的冰层温度应力。

冰的物理特性十分复杂，研究表明冰温度应力受多种因素的影响，包括冰层的热力状况、冰层的约束条件、冰层厚度等，这使得冰温度应力的估算存在诸多不确定因素。我国海冰物理力学性质[4]和动冰力[5−7]的研究相对较为系统，但针对淡水冰尤其是静冰力的研究仍然有很大缺口，国内对该问题的研究基本上采用的是现场量测冰压力的大小，经过数据分析提出相应的经验公式，将冰压力表达为冰厚和结构物几何参数的函数[8]，这些经验公式具有一定的工程应用基础，但它们的适用范围小，误差大，而且公式中一些经验参数在实际中难以确定，是一种粗略的冰力估算方法。通过建立冰层应力-应变的力学模型(本构模型)，对冰温度膨胀力的作用机理进行分析，是一类更为深入的研究思路。史庆增等[9]将层合板理论引入冰层温度应力的分析求解中，提出了一种求解对边约束冰层温度膨胀力的计算方法，并对海冰进行了计算；黄焱等[10]采用有限元方法分析了单层淡水冰的温度膨胀力，并对水工建筑物中常见的形式进行具体计算。这些基于应力-应变力学模型的冰力数值计算研究中均假定了冰的弹性本构关系。实际上，冰的本构关系与应变率明显相关，静冰温度应力发生时，弹性变形阶段只有很短的时间，因此，利用弹性本构关系求解静冰温度应力问题可能导致误差。黄文峰等[11]以红旗泡水库为例，开展不同方向上冰表面形变的原位观测，尝试改进现有的冰温度应力计算模型来预测表层应力。国外靠近北极地区的国家(瑞典、挪威、加拿大等)，对湖泊及大型水库的冰温度应力计算开展了大量的试验和理论研究，并采用了多种流变学本构模型[12−16]用于冰温度应力计算，其中大多是将两个基本的粘弹性本构模型—Maxwell模型和Kelvin模型以串联或并联的方式进行组合，然后，根据冰力原型观测数据确定模型中的待定参数，用于当地的静冰温度应力计算，是至今较常采用的静冰力计算方法。但是这些冰力模型中的参数取值在很大程度上依赖于数据获取的特定环境条件，不能直接应用于中国北方寒冷地区的静冰温度应力计算。

此外，迄今为止，国内外冰温度应力的研究多集中于湖冰和库冰，关于浅水湖塘中冰层温度应力的研究极少见于报道。与大型深水湖库相比，小型浅水湖塘的传热条件不同，从而导致冰层热力状况有所不同；另外，小型湖塘的面积更小，周边围岸对冰层的约束作用更为强烈。因此，小型浅水湖塘中冰层温度应力呈现特有的时空变化规律。

综上，国外的静冰力模型不可直接应用于中国东北地区的静冰力计算，有关塘冰温度应力的研究更是极为少见。针对上述2个研究背景，为探究我国东北地区小型静态水域的冰层温度应力变化规律，本文采用原型观测研究手段，结合理论分析及计算，以黑龙江省大庆市青花湖6号水塘为研究对象，对水塘冰层内部温度和应力进行了观测及分析，基于冰的粘弹性本构关系提出了一种塘冰温度应力计算模型，并结合实测数据回归得到了模型中的待定参数，该模型计算值与实测值吻合良好，可为寒区小型浅水湖塘及小型水库的冰温度应力计算提供参考。

1 温度应力原型观测

1.1 观测场地概况

冰情原型观测现场位于黑龙江省大庆市肇源县青花湖。青花湖，原名库里泡，位于大庆市大同区、肇州县、肇源县三地交界处，地处124.78°E～124.87°E，45.77°N～45.91°N，是安肇新河上的最终控制水体，湖区从南到北由几个相连的湖盆组成，南北长约16 km，东西最宽处约5 km，湖区水面面积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月份给湖内补水，冬季结冰期停止补水，基本没有水位变动，因此青花湖的冬季自然环境接近静态水体，提供了静冰观测的理想条件。

青花湖西端的部分水域被人为地以浆砌石挡墙划分为若干水塘进行渔业养殖，本次观测选取形状规整、塘底高程起伏不大的青花湖6号水塘为研究目标。冰温度应力主要发生在连日持续升温的条件下，根据对当地以往多年气温和冰情资料的分析，选择观测时段为2018年2月12日−3月14日，对当地气温、冰温及冰温度应力进行了持续观测。

青花湖水库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传感器薄膜区域的压力转换成电阻值的变化，从而获得压强信息，精度为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。

观测链安装前，先在冰层上沿平行于挡墙方向锯开一条贯穿的窄缝，将已经固定好温度传感器探头和压力传感器探头的观测链竖直放入窄缝之中，规定浸水后每个压力传感器的读数为该压力传感器的零点，观测链导线连接在木质支架上，木质支架固定于冰层上。冰温和冰压力观测链于2018年2月2日在青花湖水库6号水塘安装完成，当日水塘冰厚为97 cm。观测链的温度探头和压力探头随着窄缝内部冰生长被冻结到冰内，2018年2月12日窄缝冻结密实后开始进行数据采集，以保证传感器探头与冰充分接触，从而精确感知该深度处冰温及冰应力的变化。冰厚观测采用冰钻人工穿孔后用塔尺直接量测。

2 理论模型

2.1 冰的热变形

预测冰力大小的精确性很大程度上取决于冰力学模型的合理性，从冰的温度应力产生机理分析，当冰升温产生的膨胀变形受到周界约束时，冰层内部产生温度应力，温度应力的大小应恰好抵消这部分原本应该产生的温度应变，因此，冰温度应力计算主要取决于冰的热变形关系和静冰应力-应变本构关系。

由于热膨胀只产生线应变，为剪切应变为零，对于三维问题，冰层温度升高引起冰的热变形关系可表示为：

式中： ε为温度应变； ΔT/ (℃)为温度增量；α为冰的温度膨胀系数，可取α=(54+0.18T)×10−6℃−1。

如果认为膨胀系数为常数，对式(1)两边对时间求导得：

江西茶叶是以绿茶为主导，绿茶的总产量占据总产茶量的80%以上，而红茶的总产量仅占据10%左右，而其他种类的茶叶占比更是不足10%。这就在一定程度上限制了江西以绿茶出口为主的发展空间。

由式(2)可见，冰温变化率与应变速率为线性关系，因此，冰温变化引起的温度应力求解又归结为确定对应的应变速率下的冰本构关系。

2.2 流变性本构关系

静冰温度应力发生在低应变速率下(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]认为化学反应与冰的蠕变过程区别较大，这种比拟缺乏一定的解释力，并通过实测数据回归，以相对温度的指数函数形式替换Bergdahl模型中的KD，即令：

式中：T1=−1 ℃；m=1.92[13]；β均为拟合参数。

由以上推导过程可以证明，参数A即是冰弹性模量E与冰热膨胀系数 α的乘积，即：

由于A是由两个具有明确物理意义的参量组合而成，故在对式(6)拟合分析时，不再把参数A作为纯经验参数进行拟合确定，而是以某个常量代替，以避免回归分析中A、B两个参数的的相互影响。E和 α都是随冰温T变化的量，根据温度场观测结果，冰温T集中分布于−15 ℃～0 ℃区间内，取E0=4 GPa，由式(7)计算得到该温度区间内变量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商业数学软件求解，参数B由单纯形算法回归得到。

3 结果与分析

气温是影响冰层生消、冰层温度场及应力场的最主要热力因素，冰温及应力分布与气温变化紧密相关。气温观测资料(图2)表明，该观测区域于2018年2月初开始出现日平均气温回升，本研究进行数据采集的时间段为2018年2月12日−3月14日为比较稳定的连续升温阶段，符合产生冰温度应力的气温条件，观测期内冰层厚度稳定在97 cm～102 cm。下面根据实际观测数据，分别分析冰层温度和温度应力的时空分布规律。

图 2 2017年11月−2018年4月观测地日平均气温变化Fig.2 Daily mean air temperature from Nov. 2017 to Apr. 2018

3.1 冰层温度

冰层温度是冰的基本状态参量，导致温度应力产生的前提条件是温度变化与约束作用，冰温在冰层内部的变化直接影响到冰层温度应力的大小和分布。图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。通过对比应力观测值与计算值发现，尽管式(8)形式简单且只用到一个拟合参数B，但模型计算结果与观测结果吻合良好，在气温稳定回升的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节理论模型部分的介绍，可以预期参数B仅是冰的固有材料性质，然而，观测数据表明B与测点位置相关，限于本试验测点位置布置的单一性，尚无法通过观测数据确定B的成因，关于B的具体影响因素及分布规律可通过更加全面深入的观测数据进行分析。

图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日不同时刻的冰应力垂向分布廓线，由于各层冰温变化值及冰层状态不同，冰层不同深度处的温度应力也不同。冰层表面为自由面，受到的约束较小，因而产生的应力稍小，随着深度的增加冰层温度应力增大，并在冰深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，整体应力水平呈递增趋势。