水工混凝土带缝结构温度场问题的扩展有限元法

林学军，陈 安，霍中艳，陈旭东

(1.浙江省台州市港航管理局，台州 317700；2.浙江海洋大学，舟山 316000；3.郑州大学，郑州 450001)

水工混凝土带缝结构温度场问题的扩展有限元法

林学军1，陈 安1，霍中艳2，陈旭东3

(1.浙江省台州市港航管理局，台州 317700；2.浙江海洋大学，舟山 316000；3.郑州大学，郑州 450001)

长期处在温度荷载等作用下运行的水工结构或多或少都会产生裂缝，有必要对其温度场进行计算分析。鉴于扩展有限元法求解不连续问题的独特优势，文章用其求解温度场；首先分析给出了裂缝结构温度场的逼近形式，随后推导了用扩展有限元法求解不连续温度场的离散方程，其中裂缝采用等温形式进行模拟，即温度连续但其一阶导数不连续，继而建立起温度场的扩展有限元法。并以带缝十字形板的算例验证了扩展有限元法求解温度—应力问题的正确性；最后，用该方法分析了某面板堆石坝面板裂缝在冬季降温期的扩展情况，其结果与实际监测成果较一致。

水工混凝土；带缝结构；温度场；扩展有限元

从20世纪60年代起，混凝土结构在施工浇筑过程中由于混凝土水化热而产生的温度场和温度裂缝、在运营阶段由于周围环境及太阳辐射等作用产生温度荷载而导致的开裂，一直倍受世界工程界学者的关注。大量的工程实践和理论研究都表明，当温度场发生变化时，会导致混凝土体内热应力的产生，温度荷载的变化对大体积混凝土结构的应力状态有显著的影响，其导致的温度应力值可能超过其余外荷载引起的应力值，甚至比其他各种荷载产生应力的总和还要大[1]，长期处在温度荷载作用下运行的水工混凝土结构或多或少存在着裂缝，而其裂缝的存在不仅破坏了结构的连续性，还扰乱了热流的传播，从而加剧裂缝端部的应力集中现象，极易引起新的裂缝或者导致已有裂缝的扩展，对结构的安全产生不利影响，故水工混凝土结构工程中，温度场的计算分析具有非常重要的作用。

为了避免常规有限元法追踪动态裂纹时要求裂缝面与单元边界一致、网格需不断重划分等缺点，考虑到扩展有限元法对于求解不连续问题的独特优势，本文以扩展有限元法求解温度场。对于水工混凝土带缝结构而言，裂缝内介质导热系数相对周围介质较大或缝内介质温度受边界影响变化不大时，缝两侧温度相等(称等温裂缝)，可按等温裂缝处理，采用文献[2]提出的结点加强方法处理；当裂缝内介质的导热系数相对较低时(称绝热裂缝)，可通过将裂缝的导热系数取为较小值来模拟。文章首先分析裂缝结构温度场的逼近形式，然后在简述热传导问题基本原理的基础上，研究建立温度场的扩展有限元法，并以十字板算例验证了该方法的有效性，最后将其用于某面板堆石坝面板裂缝在降温作用下的扩展分析。

1 温度场的逼近形式

带缝结构温度场的XFEM逼近形式与位移场类似[2-7]，可用下式表示

式中：x为点的坐标向量；Ni(x),Nj(x)和Nk(x)为常规有限元形函数；Ti为结点温度的常规自由度，ej和fk分别为温度在各结点处的加强自由度；Φ(x)为贯穿单元结点加强函数、γ(x)是裂尖加强函数，而Φj和γk分别为二者的结点值；I为所有结点集合；J为所有与裂缝相交的单元结点集合；K为裂尖单元的结点集合。

裂缝的贯穿加强结点与位移模式中采用相同的加强函数Φ(x)，该函数连续但其一阶导数不连续。等温裂缝尖端的渐近温度场可表示为为常数，故将选为裂尖加强函数[2]。绝热裂缝尖端的渐近温度场可表示为为常数，故将选为裂尖加强函数，与扩展有限元法的位移模式原理一致[3]。

2 热传导的基本原理

2.1 热传导的方程

考虑到对于均匀、各向同性固体，由于温度升高而吸收的热量必须等于从外界流入净热量与内部的水化热之和，也即热量平衡，可知

式中：Q是单位时间单位体积内发出的热量，c和r分别为比热和密度，t代表时间，a为导温系数，a = λ / cρ。若无内热源，则温度场的热传导方程为

2.2 边值条件

为了确定所需温度场，还须知道初始条件与边界条件，二者合称为边值条件。通常认为初始瞬时温度分布为常数，即τ = 0时

边界条件一般可分为四种：

（2）第二边界条件：混凝土的表面热流量q与时间函数已知

（3）第三边界条件：假设通过混凝土表面的热流量q与其表面的温度T和气温aT的差成正比

式中：β是表面放热系数。若β→∞时，转化成第一边界条件T=Ta；若β=0时，则转化为绝热状况

（4）第四边界条件：当两种不同的固体接触时，如接触良好，则在接触面上的温度和热流均应连续，边界条件即为

2.3 XFEM离散方程

对于满足(4～5)的不考虑内热源的二维温度场问题，可根据变分原理等价于下述泛函的极值问题

将求解域划分为有限个单元，泛函Ie为单元e的求解子区域ΔR上的泛函，表达如下

温度场的计算网格及裂缝加强结点的选取与位移场和渗流场一致，加强结点仍分为裂缝贯穿加强与裂尖加强两种方式[3]，两种类型每个结点的自由度均为2个，部分双重加强的结点自由度为3个，其他常规结点的自由度为1个。

将式（9）对单元各自由度求微分，得出离散方程如下

在式(12)与(14)中，Gia= Ni，Gib= NiΨi，Gic= Niφi。

将各单元结果加以集合，则对于求解区域的全部结点可得方程组

式中：{T }为所有结点的温度阵列(包含加强温度)；[H]为整体传导矩阵，其元素Hij由与结点i有关的各单元hij与gij集合而成，[L]与[M]为时间项离散后的整体系数矩阵，其元素Lij与Mij分别由与结点i有关的各单元lij与mi集合而成，[F]为边界条件形成的温度荷载阵列；令第一类边界上各结点的温度等于给定值，求解方程(16)即可得到所有结点的温度值。

表1 工况详表Tab.1 Different conditions

图1 十字形板模型及仿真图Fig.1 Model and simulation diagrams of cross plate

3 算例分析

用第2部分介绍的方法及笔者在文献[3]中关于扩展有限元法模拟裂缝扩展的方法，分析一带有边缘倾斜裂缝的十字形板在4种不同工况下的裂缝扩展过程，各工况的温度和力学边界条件详见表1。十字形板模型的几何尺寸及边界条件如图1-a所示，其中，L = 2 m，a = 0.4 m，裂缝采用绝热加强形式，与水平轴成45°角，增量步长为0.3a；材料参数为E=2.184×105Pa，v = 0.3，α=1.67×10-5/℃。

四种工况下裂缝扩展11步后的路径及两种类型应力强度因子分别如图1-d和图2所示。为方便比较，将Ι型与Π型裂纹尖端的应力强度因子KΙ和KΠ均用K0=αE(πa0)0.5进行标准化，标准化之后的结果分别如表2以及图2所示，可发现KΙ受边界条件的影响较大； KΠ几乎不受边界条件影响，扩展一步之后趋于零；工况1下，裂纹扩展几步之后KΙ和KΠ均近似为零，说明真实条件下的裂缝应停止扩展；工况4下，裂纹扩展4步之后KΙ和KΠ趋于平稳，说明真实条件下裂纹停止扩展或处于稳定扩展状态。

图2 各工况下标准化的SIFs随裂缝扩展的变化图Fig.2 Variations of normalized SIFs by crack propagation under different conditions

图3 文献[8]边界元法计算结果Fig.3 Results of boundary element method in reference[8]

由图1-d可知，温度荷载对裂缝的扩展路径影响非常大，只加温度荷载和只加外力荷载的扩展路径夹角大于90°，温度和外力同时作用与单一外力或温度作用下的裂缝扩展路径差别也相当大，可见裂缝问题中的温度荷载是绝不能忽略的。图1-b和1-c分别给出了在工况4下裂纹扩展3步的温度场及放大1 000倍位移场。文献[8]与本文算例模型一致，采用边界元法计算，其中工况1,2,3,5分别与本文工况1,2,3,4一一对应，文献[8]各工况裂缝扩展路径图以及标准化之后的KΙ和KΠ分别如图3所示，由图可见，这与本文计算结果一致，说明了本方法的有效性。

表2 各工况下标准化的SIFs随裂缝扩展的变化表Tab.2 Variations of normalized SIFs by crack propagation under different conditions

4 工程实例

某混凝土面板堆石坝坝顶长429 m，坝顶宽度为10 m，最大设计坝高为139 m，坝顶高程2 010 m，坝体填筑分垫层、过渡料、主堆石、次堆石等区，面板共分38块，分缝长度约5 000 m，最大单块长度为218 m，其中水上面板36块；面板不等厚，厚度为0.3～0.7 m，宽度12 m，面板间设垂直缝，面板与趾板和高趾墙间设周边缝，面板混凝土设计采用525#普通硅酸盐水泥，为二级配混凝土，标号为C25；该坝的正常蓄水位为2 005 m，设计及校核洪水位分别为2 005 m及2 008.28 m。经几年运行，堆石坝混凝土面板出现较严重的裂缝问题，且绝大部分为垂直于坝轴线且纵向分布于面板上的裂缝。多次的裂缝检查发现面板混凝土受温度影响显著。面板在库水位以上部分裂缝较多，且冬季产生的裂缝较多。2011年6月最新的现场检查共发现水位以上大小裂缝135条，裂缝最长延伸至2 002.0 m高程水位以下，其中129 条裂缝为竖向裂缝，5条裂缝为横向裂缝，裂缝数量比2010年10月统计的结果增加了12条，裂缝宽度在0.02 ～0.45 mm之间。因此有必要分析温度变化对面板裂缝扩展的影响。下面以5#面板为例进行分析。

5#面板中间存在一条从2 005.5 m高程延伸至水位附近的纵向裂缝，检测发现该裂缝在2 003.5 m高程处最深，达0.11 m，宽0.15 mm。2008年初出现持续的降温过程，下面用扩展有限元法分析降温过程2003.5 m高程处裂缝的扩展深度，为此，在2 003.5 m高程处沿坝轴线的面板法向截取一个剖面，通过建立该剖面的扩展有限元法简化模型分析裂缝的扩展情况，剖面截取的示意图如图4所示，其中X向为面板宽度方向，Y向为面板厚度方向；注意到运行期面板过渡料区和垫层内的温度测值变幅较大，受外界气温影响大，而主堆石区和高程较低的过渡料区内的温度计测值的变幅则很小，故简化模型不考虑堆石区温度变化的影响。剖面的XFEM模型见图5，模型共有960个单元，1 067个结点。

由于垫层料密实，挤压墙与垫层料结合良好，运行期的挤压墙压碎后性能已经基本与垫层相同，因此，在计算温度应力时仅考虑挤压墙与面板间的接触问题，挤压墙与垫层之间的接触问题不再考虑。通过在面板与挤压墙之间设置无厚度Goodman接触单元模拟二者的接触问题。

图4 面板法向剖面的截取图Fig.4 Profile normal of the concrete panel

图5 面板剖面的网格模型(单位: m)Fig.5 Mesh model of the panel profile

据面板坝所处地区的气候特点，由气温观测资料，选取旬平均气温最低的2008年1月下旬为典型工况之一进行分析，气温变化过程如图6所示；由于相邻面板之间温度场影响较小，计算温度场时，面板左右侧取为绝热，即仅考虑热流沿面板法向流动，不考虑其纵向流动；面板混凝土的材料参数见表3，混凝土与空气的热交换系数β=2 009.28 kJ/m.d.℃，断裂能Gf= 120 N/m。

接触面单元劲度系数的选取[9]：（1）法向劲度系数，当接触面受压时，为了防止两边二维单元在接触处重叠嵌入，应取较大数值，本节取107 kN/m3；若算出法向应力为拉应力，假定接触面不能受拉，则取较小值，本节取10 kN/m3，使算出的拉应力可忽略不计（2）切向劲度系数用适合增量法的经验公式，每级荷载更具接触面单元的应力状态均用公式计算。

图6 典型工况下气温过程线Fig.6 Temperature changes in typical working conditions

表3 面板混凝土材料参数表Tab.3 Material parameters of face concrete

通过对整个工况的温度场计算及温度计监测资料可知，面板沿厚度方向从表面至底部受气温影响逐渐变弱，面板表面温度与气温相同，最大温差为-5.5℃，出现在2008年1月30日。

采用增量法求解冬季降温期工况下(如图6所示)面板裂缝的扩展过程，初始温度场取年平均气温8.5℃，认为此时面板内部无应力作用，所谓的增量是指某时刻温度对比初始温度场的增量，正负均可；若温度增量为正，表示温升，此时温度应力相比前一日有所减小（拉应力为正），反之则温度应力增加。按照文献[3]进行迭代计算，以时刻tn的温度场为初始条件，得出时刻tn+1=tn+Δt的温度场，以此类推，不断求解新的温度场和应力场，并以最大主应力准则判断裂缝在各时间步是否扩展，最终得出整个降温期的扩展路径。时间步取为1 d；通过反复试算，裂缝增量长取为0.01 m；计算得出裂缝的扩展路径如图7所示；2008年1月29日，裂缝达到最大深度0.27 mm；图8为2008年1月29日裂缝稳定后裂尖附近X向的局部应力放大图。

由图7和图8可知，温降拉应力使得面板缝向垫层延伸扩展，拉应力在裂尖附近有一定的应力集中现象，随着温度降低，裂缝不断扩展，裂缝扩展深度为0.27 m，比检测值大0.03 m(钻芯取样的检测值为0.24 m)，可能是由于没考虑三维情况及面板内的钢筋的缘故，扩展路径与监测结果较一致。

图7 裂缝扩展路径Fig.7 Crack path

图8 裂缝稳定后(2008年1月29日)裂尖附近的X向应力图(单位: MPa)Fig.8 Distribution of stress near the crack-tip in X direction when crack being stable on January 29, 2008

5 结论

鉴于温度荷载对水工结构此类大体积混凝土结构的影响较大，易导致裂缝进一步扩展，产生不利影响，考虑到扩展有限元法处理不连续问题的独特优势，故本文采用扩展有限元法求解水工带缝结构的温度场，推导了用扩展有限元法求解不连续温度场的离散方程，其中裂缝采用等温形式进行模拟，即温度连续但其一阶导数不连续，并以带缝十字形板的算例验证了扩展有限元法求解温度-应力问题的正确性和有效性，最后将该方法用于某面板堆石坝面板裂缝在冬季连续降温作用下的裂缝扩展分析，其裂缝扩展深度、路径与实际监测结果较一致。

[1] 龚召熊,张锡祥,肖汉江，等. 水工混凝土的温控与防裂[M]. 北京:中国水利水电出版社,1999.

[2] Sukumar N, Chopp D L, Moes N, et al. Modeling holes and inclusions by level sets in the extended finite-element method[J]. Computer Methods in Applied Mechanics and Engineering, 2001,190(46-47): 6 183-6 200

[3] 霍中艳,郑东健.基于扩展有限元法的粘聚性裂缝模型的混凝土梁断裂过程模拟[J]. 计算机辅助工程,2010,19（4）:29-33. HUO Z Y, ZHENG D J.Simulation of concrete cracking using cohesive crack model based on XFEM[J]. Computer Aided Engineering,2010,19(4):29-33.

[4] Dolbow J, Moes N, Belytschko T. An extended finite element method for modeling crack growth with frictional contact[J]. Computer Methods in Applied Mechanics and Engineering, 2001, 190(51-52): 6 825-6 846.

[5] Khoei A R, Nikbakht M. An enriched finite element algorithm for numerical computation of contact friction problems[J]. International Journal of Mechanical Sciences,2007, 49(2): 183-199.

[6] Moes N, Cloirec M, Cartraud P, et al. A computational approach to handle complex microstructure geometries[J]. Computer Methods in Applied Mechanics and Engineering, 2003, 192(28-30): 3 163-3 177.

[7] Sukumar N, Chopp D L，Moes N, et al. Modeling holes and inclusions by level sets in the extended finite-element method[J]. Computer Methods in Applied Mechanics and Engineering, 2001, 190(46-47): 6 183-6 200.

[8] Prasad N N V, Aliabadi M H. Incremental crack growth in thermoelastic problems[J].International Journal of Fracture, 1994, 66:45-50.

[9] 殷宗泽. 土工原理[M]. 北京:中国水利水电出版社,2007.

Temperature field analysis of harbor engineering structure with cracks based on XFEM

LIN Xue-jun1, CHEN An1, HUO Zhong-yan2, CHEN Xu-dong3
(1. Taizhou Port and Shipping Management Bureau,Taizhou 317700,China; 2. Zhejiang Ocean University,Zhoushan 316000,China; 3. Zhengzhou University,Zhengzhou 450001,China)

For the long-term exposure to temperature load, certain parts of a concrete structure have more or less cracks during its construction and operation period. It is necessary to calculate and analyze the temperature fi eld. For the diffi culty of applying classical fracture criteria to the actual hydraulic engineering and simulating the process of cracking by conditional FEM, the XFEM was introduced into the analysis of the temperature fi eld in hydraulic structures in this paper. Firstly, the approximation form of the temperature fi eld of cracked structure was analyzed and provided, and then the discrete equation for solving discontinuous temperature fi eld was deduced, where the crack was simulated by isothermal form, namely, the temperature field is continuous but its first derivative is not continuous. The temperature fi eld by XFEM was established. And a numerical example of a cross shaped plate with slit was used to verify the validity of proposed method. Finally, this method was used to analyze the expansion of the cracks in the face slab of a concrete face rock-fi ll dam in winter. The results are consistent with the actual monitoring results.

harbor engineering concrete; structure with crack; temperature fi eld; XFEM

TU755；O 241.82

A

1005-8443(2017)03-0274-07

2016-11-25；

2017-02-16

浙江省交通厅科技计划项目：港口码头健康检测评估、修复加固技术水文水资源与水利工程科学国家重点实验室开放基金（2015491411）

林学军（1966-），男，浙江台州人，高级工程师，主要从事港口航道工程建设与管理方面工作。

Biography:LIN Xue-jun(1966-),male,senior engineer.