基于有限单元法寒区渡槽静冰压力分布特征研究

李可也，袁 理

（中水淮河规划设计研究有限公司，合肥 230601）

高海拔高寒地区冬季气温较低， 且白天太阳辐射强烈，近地面的气温会迅速上升，日温差较大。由于高原低温期持续时间较长，水库、河流、渠道内极易结冰。研究表明，冰在气温升高时体积膨胀［1］，会对约束其膨胀的水工建筑物产生温度膨胀力， 又称静冰压力。 较大的静冰压力会对水工建筑物产生严重危害， 为研究高寒地区较大日温差下的静冰压力分布规律， 本文选取冰盖形状和约束形式都比较简单的渡槽进行仿真分析， 通过对其温度场与应力场的耦合计算，研究渡槽内冰盖在不同温升率、不同冰厚、不同约束形式下的静冰压力分布。

1 静冰压力计算原理

静冰压力主要与冰盖的温度分布、温升率、冰盖厚度等因素有关， 运用温度加载法对静冰压力进行仿真研究时，由于加载速率较慢，因此冰表现出一定的黏弹性［2］，黏弹性材料的温度变化与材料变形的关系可表示为：

式中ε 为冰的线应变；α 为冰的温度膨胀系数，一般认为淡水冰的温度膨胀系数为常量；Δθ 是温度的增量。

对式（1）的时间分量求导，得：

由式（2）可知，冰的应变速率与温升率呈线性关系，静冰压力的大小受温升率的控制。

静冰压力求解过程中首先需要计算出冰盖的温度场， 冰盖内部温度场的变化主要由冰表面的大气温度变化引起的，在有限单元法中，这是一个瞬态传热问题， 根据能量守恒原理， 瞬态热平衡可表达为（以矩阵形式表示）［3］：

式中[K]为热传导矩阵，包含导热系数、对流系数及辐射率和形状系数；[C]为比热矩阵，考虑系统内能的增加；{T}是节点温度向量；{T·}为温度对时间的导数；{Q}为冰点热流向量，包含热生率。

计算渡槽内冰盖温度场时， 认为淡水冰盖下表层与水的接触面的温度为其融化点0 ℃［4］，冰盖上表层温度按第三类边界条件确定， 即冰层上表面的温度与大气温度相同。第三类边界条件可表示为：

式中 Tf为热流介质温度；Γ 为物体边界；ζ 为热换系数。

求得冰盖温度场后， 将有限元模型中各节点的温度作为荷载引入应力场计算中， 通过定义材料的热膨胀系数α 将温度场与应力场进行耦合计算，最后得到冰盖内的应力分布规律， 作用于边界上的冰的温度膨胀力即为冰对水工建筑物的静冰压力［4］。

通过原型观测，计算冰盖应力场时，做以下两点假定［5］：

（1）计算初始时，冰盖内部无温度膨胀产生的应力，边界上的静冰压力为零。

（2）认为冰盖边界与渡槽冻结，仿真计算时，约束类型为四周全约束。

2 直立墙渡槽内的静冰压力

为研究高寒地区静冰压力分布规律， 本文以地处高原地区平均海拔2300 m 的某县某直立墙渡槽为例进行了静冰压力的仿真研究，据气象资料查得，该县近年来测得的年最低温度为-20 ℃，最大日温升达2 h 升高10 ℃。冰盖各项参数取值如下：

表1 冰盖计算参数

2.1 冰盖温度场分布

取直立墙渡槽内2 m×8 m×0.4 m 的冰盖为研究对象，为取得较精确的结果，将冰盖沿厚度分为20层，模型共划分单元32000 个，节点139041 个，冰盖有限元模型如图1。首先计算冰盖温度场分布，设冰盖上表面度为-20 ℃，且在2 h 升高10 ℃，冰盖下表面温度为0 ℃，计算初始时，认为冰盖已形成稳定的温度场。将稳态温度场作为初始条件进行迭代，最终计算得到温升2 h 过程的瞬态温度场分布， 稳态温度场云图和沿程瞬态温度场曲线如图2、图3。

图1 冰盖三维有限元计算模型

图2 冰盖稳态温度场分布云图

图3 冰盖沿程瞬态温度场分布曲线

由图2 和图3 可看出，渡槽内的稳态温度场呈线性分布，随着时间的增加，冰盖上表面温度的升高，冰盖温度场产生非线性变化，但由于冰的导热性较差，温差的升高集中在冰盖上部1/3以上处，且从表层向下变化程度逐渐减小，瞬态温度场在约15 cm 处以下基本与稳态温度场重合。

2.2 静冰压力分布

得到冰盖的温度场之后， 将热单元转换为结构单元。由于假定冰层内部初始应力为零，这就要求给材料赋予线膨胀系数，设定材料参考温度时，将稳态温度场计算时求出的单元温度值分布赋予应力场分析时对应的单元。 然后将瞬态温度场计算得到的各步温度作为荷载施加到应力计算的模型中， 经过迭代得出各节点的应力分量［6］。由于所求结果为压应力分布，将Y 方向的应力作为输出结果进行分析，直立墙渡槽内0.4 m 厚冰盖的静冰压力云图分布如图4～图7，冰盖沿程瞬态应力场分布曲线如图8。

图4 0.5h 时冰盖的应力场分布云图

图5 1h 时冰盖的应力场分布云图

图6 1.5h 时冰盖的应力场分布云图

图7 2h 时冰盖的应力场分布云图

图8 冰盖沿程瞬态应力场分布曲线

由图8 可看出，随着冰盖表层温度的缓慢上升，冰盖产生的温度膨胀力受到周围约束的作用产生静冰压力。静冰压力沿程分布不均，集中在冰盖上部1/3以上处，下部趋近于零，表层静冰压力最大，静冰压力随着温度的不断升高持续增大，2 h 的最大静冰压力值为268.7 kPa，这与温度场仿真的结果相对应，即温度变化最大处的静冰压力也最大。

2.3 温升率对静冰压力的影响

为了研究温升率对静冰压力的影响，假定大气温度在 2 h 内从-20 ℃升高至-15 ℃，其他条件不变， 得到冰盖的应力场仿真结果，2 h时两种温升率下的Y 方向静冰压力沿程分布曲线如图9。

图9 两种温升率下2h 时冰盖沿程静冰压力分布曲线

由图9 可看出，两种温升率下的静冰压力均集中于冰盖1/3 以上处，静冰压力最大处同样位于冰盖表层。温升率的大小决定了静冰压力的大小，随着温升率的增大，静冰压力也变大，符合之前的理论分析。

2.4 冰厚对静冰压力的影响

为研究不同冰盖厚度对静冰压力的影响， 根据该县当地真实的冰清资料进行假定， 当冰盖厚度分别为0.3，0.4，0.5，0.6，0.7，0.8 m， 按2 h 温升10 ℃，其他条件不变时，冰盖表层最大静冰压力值分布如图10。

图10 不同厚度冰盖表层静冰压力分布曲线

由图10 可看出，冰盖表层静冰压力的最大值随着冰盖厚度的增大而增大。且由仿真分析可知，不同厚度的冰盖，其静冰压力仍集中作用于冰盖表层值冰厚15cm处，说明冰盖表层的温升产生的静冰压力仅对冰盖上部产生比较大的影响，对冰盖下部影响很小。

3 结语

（1）采用有限元法对冰盖的静冰压力进行仿真计算具有相当精度， 在进行高寒地区水工建筑物设计时， 其表层最大静冰压力仿真结果可作为估算冰盖静冰压力极值的参考。

（2）冰盖瞬态温度场变化集中于冰盖上部15cm以上处。

（3）温升率对冰盖表面静冰压力有较大影响，静冰压力随温升率增加而增加。

（4）冰盖厚度对表层最大静冰压力值有一定影响，但静冰压力仍集中于冰盖上部15cm 以上处。