考虑变温效应的重力坝裂纹尖端场有限元分析

Finite Element Analysis of Crack-Tip Fields for a Cracked Gravity Dam with the Effect of Temperature Change

摘要: 本文基于ABAQUS对刘家峡水电站含有裂纹的非溢流坝段进行有限元模拟分析，考虑由于季节变化引起的变温效应，获得了裂纹诱导的热应力场。结果表明，裂纹前沿应力场随着温度的变化而变化，与不考虑温度场的情形比较，变温引起的应力集中能够导致裂纹扩展，所得结论能为大坝的寿命估计提供一定的参考依据。

Abstract: Based on ABAQUS, finite element analysis is carried out for the non-spillway section with cracks in Liujiaxia hydropower station. And the thermal stress field induced by cracks has been obtained by considering the effect of temperature change with the seasons. The obtained results show that with the change of temperature, the crack-tip fields are changed. As compared with the case without considering the temperature filed, the stress fields induced by temperature change could lead to crack growth. The observations provide a certain basis for evaluating the service life of the gravity dam.

1. 引言

由于混凝土重力坝坝体材料和大体积特点，施工中受温度影响等因素，坝体表面不可避免会产生微裂纹，特别是一些纵向裂纹 [1] [2]。这些初始微裂纹的存在会在一定程度上降低大坝的安全系数。在运行过程中，随着库水位的升高及外部温度场的变化，坝踵区拉应力值增加，微裂纹将会扩展并造成坝踵开裂，影响大坝的安全性。根据近些年来的研究显示，影响重力坝应力的因素除了自重、水压力、地震等之外，温度应力也是重要组成部分之一。例如：在三峡大坝孔口附近位置，温度应力的影响作用程度大于坝体自重和内外水压力之和 [3]。温度会对坝体造成影响，主要是由于外界温度与坝体内部温度存在温度差。温度梯度的存在使混凝土受力不均，一旦超过了混凝土的最大抗拉强度则会引起初始裂纹的继续扩展 [4]。

近年来，国内外许多学者都对混凝土温度场进行了不少研究分析。戴跃华等 [5] 通过模拟分析隔河沿大坝运行期内月平均日平均变化的温度场，得到了在温度场作用下坝体位移应力的变化规律，并提出尽量避免冬季温度持续下降时水库在低水位运行，夏季长时间高温环境下在高水位运行的建议；袁自立等 [6] 通过模拟石漫滩水库运行期的温度场和温度应力场，从而发现了引起坝体横向裂纹产生和扩展的主要因素及其危害；张子明等 [7] 利用裂纹带模型，对大体积混凝土的温度应力和温度开裂进行了计算；沈健 [8] 对某斜拉桥塔即有裂纹在温度疲劳荷载下的扩展概率进行分析，确立带预裂纹混凝土疲劳应力与疲劳寿命的对应关系；唐世斌等 [9] 通过模拟在升温情况下混凝土裂纹的产生及扩展过程，揭示了主裂纹形成的根本原因是裂纹扩展所伴随的局部化效应；郭生根 [10] 模拟分析了混凝土受热膨胀至产生裂纹的过程，获得了混凝土最大主应力与外界温度变化的函数关系；刘杏红等 [11] 利用伽辽金无网格法编制了分析程序，模拟了混凝土块在温度应力作用下的裂纹的产生和扩展过程，得到了在均匀温降和基础强约束条件下大体积混凝土试块裂纹的分布规律；王亮明 [12] 采用ANSYS，模拟分析寒潮期间刘家峡水电站非溢流坝段坝体温度场和温度应力变化情况，得出冬季温度所引起的拉应力十分接近混凝土抗拉强度的结论。因此研究在正常运行条件下，变化的温度场对重力坝表面裂纹的影响，对评估混凝土重力坝的安全超载能力具有重要意义。

本文通过ABAQUS平台，对刘家峡水电站含有裂纹的情形进行有限元计算，分析了季节变化导致的温度场对坝体裂纹尖端场的影响，所得结果能够为大坝的寿命估计提供一定的理论基础。

2. 热传导方程

基于Fourier热传导理论，热流q与温度梯度 $\partial T/\partial x$ 成正比，但热流方向与温度梯度方向相反，即

$q_x=-\lambda \frac{\partial T}{\partial x}$ (1)

其中 $\lambda$ 为导热系数( $kJ/{\text{m}}^{\text{2}}\text{s}℃$ )。由热量的平衡原理，温度升高所需吸收的热量必等于从外界流入的净热量与内部产生热量的总和。化简后的二维热传导方程如下：

$\frac{\partial T}{\partial t}=\alpha \left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}\right)+\frac{\partial \varphi }{\partial t}$ (2)

其中 $\alpha$ 为导温系数 $\left({\text{m}}^2/\text{s}\right)$ ；t为时间；φ为混凝土材料的绝热温升(℃)。在初始条件下，即 $t=0$ 时， $T\left(x,y,0\right)=T\left(x,y\right)$ ，其中 $T_0\left(x,y\right)$ 为已知的初始温度。

本次模型所采用的边界条件如下：

第一类边界条件： $T\left(t\right)=f(\; t\; )$

其中 $f\left(t\right)$ 为定义的边界条件。此条件可用于坝体表面与库水接触的部分，此时坝体表面混凝土的温度与所接触的水温相等。

第三类边界条件： $-\lambda \frac{\partial T}{\partial x}=\beta \left(T-T_a\right)$

式中β为表面放热系数( $kJ/{\text{m}}^{\text{2}}\text{s}℃$ )，当 $\beta =0$ 时，即为绝热条件。此条件可用于坝体表面和大气接触的部分，此时坝体表面的热流量与坝体表面混凝土温度T和大气温度T_a之差成正比 [13]。

3. 计算模型及参数

3.1. 计算模型

刘家峡水电站简化后的坝段模型几何尺寸如下图所示。以坝踵为原点、顺水流方向为X轴正向、垂直方向为Y轴建立坐标系。其中坝高H为79 m，坝顶宽度为16 m，坝底宽57.6 m，下游坝坡1:0.65，从坝底到坝高64 m处 [12]。坝基计算范围为往上游、下游均取近1倍坝高即100 m，基岩计算深度也同样取为100 m。具体断面尺寸及初始裂纹所在位置如图1所示，坝体材料性能参数见表1。

Figure 1. Non-spillway section model (unit: m)

图1. 非溢流坝段断面尺寸(单位：m)

Table 1. The table of material performance parameters

表1. 材料性能参数表

3.2. 计算工况

刘家峡水电站在一月份时最大库水位为75 m，在七月份时最大库水位为50 m [12]。为了便于比较变化的温度场对坝体初始裂纹张开宽度、尖端应力等影响，以下分三种工况进行对比：1) 不施加温度边界条件，对含有初始裂纹的坝段进行应力分析；2) 以1月份寒潮为温度边界条件，对含有初始裂纹的坝段进行应力分析；3) 以7月份气温日变化为温度边界条件，对含有初始裂纹的坝段进行应力分析。此次分析所设置的初始裂纹均为开口0.01 m、长5 m的水平裂纹，示意图如图2(a)所示。裂纹分别处于迎水面上坝底往上20 m、40 m、60 m的位置，依次编号为L20、L40、L60。同时，在文献 [12] 的基础上，对工况二、三的初始裂纹附近的应力场进行比较，研究在相同温度场下，初始裂纹对坝体的影响。

如图2(b)为网格划分的示意图。对于工况一，采用CPS4型网格：四结点双线性平面应力四边形单元；对于工况二、三，采用CPS4T型网格：双向线性位移，双向线性温度，四结点平面应力热耦合四边形单元。

(a) 初始裂纹断面示意图(单位：m) (b) 网格划分示意图

Figure 2. Sketch diagram of model

图2. 计算模型示意图

3.3. 温度条件

刘家峡水电站位于我国西北高原区域，属于大陆性气候，寒潮多发生在1月份，日气温在48~72小时内急剧下降。通过用正弦函数，设在t₁时遇到寒潮，48小时内的气温变化可用下式计算 [13]：

$T_a=T_m-A_{48}\sin \left[\frac{\pi \left(t-t_1\right)}{2Q}\right],t\ge t_1$ (3)

式中：T_m为月平均气温，在1月份取−4.9℃；A₄₈为寒潮中的气温变幅，此处取10℃；Q为寒潮中降温历时，此处为48小时；t为温度变化时间；t₁为寒潮起始时间。

7月份为当地年最高气温月，昼夜温差较大。通过用余弦函数，日气温变化可用下式计算 [13]：

$T_a=T_m-A_{24}\cos \left[\frac{\pi }{12}\left(t-t_2\right)\right],t\ge t_2$ (4)

式中：T_m为月平均气温，在7月份取22.4℃；A₂₄为24小时内气温变幅，此处取17℃；t为温度变化时间；t₂为气温最高时间。

在坝体与库水接触部位，坝体混凝土表面温度可采用第一类边界条件进行计算。综合库水位、气象资料和水温分布特点分析，任意深度的水温计算公式为 [14]：

$T_m\left(y\right)=c+\left(T_s-c\right){\text{e}}^{-0.04y}$ (5)

$c=\frac{T_b-T_{s}g}{1-g}$ (6)

$g={\text{e}}^{-0.04H}$ (7)

$T_s=T_a+\Delta b$ (8)

式中：T_b为库底年平均水温，按照朱伯芳院士提出的方法，处于寒冷地区且深度在50m以上的水库T_b一般取6℃；T_s为水库表面水温；T_a为气温；Δb为温度增量，主要受到日照的影响，冬天取1℃，夏天取3℃；H为水库水深。

在坝体与空气接触部位，坝体混凝土表面温度可采用第三类边界条件。

4. 有限元计算成果及结论

4.1. 温度场分析

在1月份48小时寒潮期间，气温持续下降10℃。图3表示寒潮期间初始裂纹L20模型的坝体温度场分布情况，其中(a)、(b)各代表寒潮历时1个小时和48小时后的情况。由图可知，随着温度的下降，与空气直接接触部分的气温也在持续下降，在48小时内温度由−4.9℃下降到−14.9℃。温度大幅度变化的范围仅限于混凝土表层部分，对坝体内部的影响相比较小。由于裂纹的存在，增加了坝体混凝土表层面积，使坝体内部温度场受外界气温变化影响的范围增大。其中，裂纹位置越靠近坝底，影响坝体内部温度场不均匀分布越明显。

在7月份24小时气温日变化期间，气温变化幅值为17℃。图4表示气温日变化期间初始裂纹L20模型的坝体温度场分布情况，其中(a)、(b)各代表24小时内气温最低和最高时刻的温度场情况。由图可知：随着温度的变化，与空气直接接触部分的气温也在呈现余弦式变化，在凌晨2:00点时气温最低，中午14:00点气温最高，在24小时内温度由30.07℃下降到14.70℃，而后上升回30.07℃。与气温直接接触的坝体混凝土表面部分，受气温影响产生的温度变化幅值与坝体内部相比更大些。由于上游水深为50 m，L60裂纹此时已暴露在空气中，直接受到气温的影响。对于位于水下的初始裂纹，裂纹位置越靠近坝底，影响坝体内部温度场分布越明显。

(a) 寒潮1个小时后 (b) 寒潮48个小时后

Figure 3. Nephogram of cold wave temperature field in Jan

图3. 1月份寒潮温度场分布云图(L20)

(a) 气温最低 (b) 气温最高

Figure 4. Nephogram of daily temperature field in July

图4. 7月份日温度场分布云图(L20)

4.2. 裂纹尖端应力场

下面主要分析裂纹尖端场。

在1月份寒潮工况下，部分裂纹模型的应力云图见图5、图6，其中(a)、(b)各代表寒潮历时1个小时和48小时后的情况。

在X方向，裂纹附近部位的压应力由裂纹开口处朝尖端逐渐增大。对比不同的裂纹位置情况可得，裂纹位置高程越高，越靠近水面，附近的X方向压应力越小。随着温度的降低，压应力在不断减小，压应力最大值会出现在裂纹尖端位置，最小值由坝体内转移到裂纹中端。在L60裂纹模型中，裂纹中端434号结点X方向的应力由−0.422 MPa增加到0.14 MPa，这说明出现了该节点附近的应力状态随着温度的降低，由压应力转变为拉应力。在以后坝体工作的过程中，若遇到更强降温的不利工况，该部位的拉应力一旦大于混凝土抗拉强度，就会引起裂纹的扩展，并向坝体内部发展。而与工况一相比，当不考虑温度场时，裂纹附近均处于受压状态。故对含有初始裂纹的坝体分析时，温度场也是一个重要考虑因素。

(a) 寒潮1个小时后 (b) 寒潮48个小时后

Figure 5. Nephogram of X-direction stress in Jan (L60)

图5. 1月份X方向应力分布云图(L60)

在Y方向，同样，裂纹附近部位的压应力由裂纹开口处朝尖端逐渐增大。对比不同的裂纹位置情况可得，裂纹位置高程越高，越靠近水面，附近的Y方向压应力越小。随着温度的降低，压应力也在不断减小，压应力最小值会出现在上游坝面裂纹上方附近，最大值在裂纹尖端位置。

与工况一相比，当不考虑温度场时，裂纹附近均处于受压状态，且随着裂纹开口处朝尖端逐渐增大，Y方向的最小压应力为0.505 MPa，最大压应力为1.235 MPa。在温度场的影响下，裂纹附近的压应力逐渐减小。随着温度降低，最大压应力由1.722 MPa减小到1.606 MPa，但仍比不考虑温度场的模型大。

(a) 寒潮1个小时后 (b) 寒潮48个小时后

Figure 6. Nephogram of Y-direction stress in Jan (L20)

图6. 1月份Y方向应力分布云图(L20)

对于在7月份日气温变化工况下，部分裂纹模型的应力云图见图7、图8，其中(a)、(b)各代表24小时内气温最低和最高时刻的情况。由于7月份上游水深为50 m，此时L60裂纹模型中的裂纹已暴露在空气中。

在X方向，裂纹附近部位的压应力由裂纹开口处朝尖端逐渐增大。对比不同的裂纹位置情况可得，裂纹位置高程越高，越靠近水面，附近的X方向压应力越小，在L40和L60裂纹模型中已经出现了拉应

(a) 气温最低 (b) 气温最高

Figure 7. Nephogram of X-direction stress in July (L40)

图7. 7月份X方向应力分布云图(L40)

力。随着温度的降低，各模型应力数值都在不断变化，压应力最大值均出现在裂纹尖端位置。对于L20裂纹模型，压应力最小值在坝体内部；对于L40和L60裂纹模型，X方向的最大应力均为正值，这说明出现了该节点附近的应力状态为拉应力。在温度场的影响下，对于L40裂纹模型，裂纹附近X方向的最大应力出现在坝体内部，由0.068 MPa增加到0.073 MPa；最小应力为裂纹尖端处，由−1.308 MPa增加到−1.332 MPa。

与不考虑温度场工况相比，裂纹附近X方向的应力状态，随着裂纹开口处朝尖端逐步减小，最大应力为0.071 MPa，最小应力为−1.291 MPa。这说明，在坝体工作的过程中相应位置X方向的应力会随着温度的升高而减小，降低而增加，并且在拉、压应力状态之间来回切换。

(a) 气温最低 (b) 气温最高

Figure 8. Nephogram of Y-direction stress in July (L60)

图8. 7月份Y方向应力分布云图(L60)

在Y方向，同样，裂纹附近部位的压应力由裂纹开口处朝尖端逐渐增大。对比不同的裂纹位置情况可得，裂纹位置高程越高，越靠近水面，附近的Y方向压应力越小，在L60裂纹模型中已经出现了拉应力。随着温度的降低，各模型应力数值都在不断变化。其中，压应力最大值始终处于裂纹尖端位置。对于L20和L40裂纹模型，压应力最小值在上游坝面裂纹开口处；对于L60裂纹模型，X方向的最大应力为正值，这说明出现了该节点附近的应力状态为拉应力。

与不考虑温度场工况相比，对于L60裂纹模型，当不考虑温度场时，裂纹附近Y方向的应力状态，随着裂纹开口处朝尖端逐渐减小，最大应力为0.098 MPa，最小应力为−1.648 MPa。在温度场的影响下，裂纹附近Y方向的最大应力出现在上游坝面裂纹上方，由0.098 MPa增加到0.099 MPa；最小应力为裂纹尖端处，由−1.645 MPa增加到−1.640 MPa。这说明，在坝体工作的过程中相应位置Y方向的应力会随着温度的升高而减小，降低而增加。

5. 结果与讨论

从以上分析可以看出，在1月份寒潮工况下，对于含初始裂纹的模型，X方向的压应力最大值在裂纹尖端处，Y方向的压应力最小值在上游坝面裂纹上方，都随着温度降低而减小。在7月份日气温变化工况下，对于含初始裂纹的模型，X方向的压应力最大值在裂纹尖端处，在Y方向的应力最大值在上游坝面裂纹上方，均随着温度降低而增加。

对于不考虑初始裂纹的模型，各方向上的应力情况受温度影响较弱，相同位置的压应力随温度变化而变动的幅度较小。且在相同时期，X方向上的压应力要比含初始裂纹的模型小，Y方向上的压应力要比含初始裂纹的模型大。

结合两种工况分析，对于处在水面波动区域的裂纹尖端，该处的应力会随着温度的变化而在拉压状态下切换。在实际的工程运行中，遇上更不利的变温、较大的水位波动情况仍会存在。因此，在后期坝体工作的过程中应采取相应措施，尽量避免裂纹扩展，降低坝体的整体性，以免危害到大坝的安全。

基金项目

国家自然科学基金资助(No. 11872155)。