基于应力强度因子的混凝土大坝诱导缝变形数值模拟

顾 双，殷宪柱，程 炀，鲍海兵

(1.江苏省骆运水利工程管理处，江苏 宿迁 223800；2.宿迁市水务勘测设计研究有限公司，江苏 宿迁 223800；3.扬州市勘测设计研究院有限公司，江苏 扬州 225500；4.江苏省水利建设工程有限公司，江苏 扬州 225000)

20世纪70年代，美国首次提出混凝土筑坝技术，后于1981年，在日本建造了第一座重力坝——岛地川坝。建立的混凝土重力坝，水泥使用量较小、混合掺量材料量较高，降低了工程成本的同时，保证了大坝的抗冲击性能。我国不甘落于人后，在20世纪80年代，展开了对混凝土大坝筑坝技术的研究，在1986年，建造了国内第一个混凝土重力坝。就目前的技术而言，为了进一步加强大坝的安全性，在扩大混凝土大坝建设范围的基础上，提出通过分析大坝诱导缝受力变形状态，加强大坝的使用安全[1]。一些学者通过分析诱导缝布设形式对拱坝应力变位的影响，计算坝体拉应力与压应力最大值、坝体位移等参数，实现对混凝土大坝诱导缝的模拟分析[2]。还有学者应用大型有限元分析软件ANSYS，模拟混凝土大坝诱导缝，通过三维有限元应力仿真分析法，分析诱导缝开裂情况，实现对混凝土大坝的可行性研究[2- 3]。经过实际应用测试发现，这些方法的使用效果并不理想，因此需研究全新的数值模拟方法，基于应力强度因子分析混凝土大坝诱导缝的受力变形情况。该方法为模拟方法，通过模拟混凝土大坝中所有有效数据，实现对应力值、受力值等参数的分析。

1 数值模拟方法

1.1 计算诱导缝应力强度因子

应力强度因子是一个有效状态参量，表示裂纹尖端前缘应力场强弱程度，该数据不只是表示应力应变数值的变化情况，还会表示整个场的能量值，是构建诱导缝等效强度双参数受力变形模型用到的必要数值。根据弹性断裂力学原理可知，应力强度因子与坐标无关，只是反映应力场奇异性的强弱，其量纲为力和长度的乘积[4]。应力强度因子Qi对于不同开裂类型的大坝，具有不同的通式：

(1)

式中，QI—张开型大坝；QII—滑动型大坝；QIII—撕开型大坝；θ、β、γ—对应大坝裂纹的几何形状因子；φ—拉力；η1、η2—内切应力和外切应力。

在传统数值模拟方法中，进行有限元计算时，应力强度因子影响模型模拟大坝诱导缝的受力变形过程，因此本次研究的数值分析方法，计算诱导缝应力强度因子时，采用复变函数修正法，获取裂纹尖端应力强度因子的具体数值，计算公式为：

(2)

式中，n—其他类型大坝；N—大坝类型；μ0—裂尖坐标，通常情况下该值的取值为1或-1[5- 6]；ε(*)—复变函数的分离结果。

上述公式是面对无限板单孔边裂纹结构的计算公式，下列公式应用于无限板，多孔边裂纹结构的应力函数：

(3)

式中，εj—孔口j的复变应力函数；μj—与第j个函数对应的坐标[7]。则应力强度因子的计算公式为：

(4)

根据需要模拟的大坝实际情况，选择对应的计算公式，获取诱导缝应力强度因子。

估摸时间差不多了，玉敏收起搓澡巾，说姑妈，好了。姑妈说我再帮你搓会？玉敏赶紧摆手，说你要折杀我啊，你是姑妈，我哪能劳驾你！再说我的背不用搓，每次都是小虫搓的。说完不好意思地笑了。两人又去冲了淋浴，半小时后回到休息室，闲聊了会。姑妈问起杨律师的事，玉敏说杨律师和那客户协商未成，已向法院提起诉讼了，过几天就开庭。姑妈说这官司有把握么，玉敏说杨律师说有，但没把握能拿到钱。姑妈叹息。玉敏又说杨律师说了，如果那客户有财产，可以执行的。姑妈不懂法律，哦了一声，说那就催杨律师抓紧吧。

1.2 构建诱导缝等效强度双参数受力变形模型

根据传统数值模拟方法的实验研究结果可知，混凝土大坝受力变形过程可分为：初始起裂、裂缝扩展以及破坏失稳三个不同的阶段。

(5)

式中，t0—微裂缝生成区段。

根据荷载和初始裂缝长度，计算应力强度因子Qi的起裂韧度，当荷载量达到最大值时，应力强度因子达到应力临界值Qmax，此时的裂缝会破坏失稳，断裂区长度λ达到临界值Δλ[9]。图1为裂缝区应力直线分布的断裂区长度、有效裂纹长度示意图。

图1 断裂区长度、有效裂纹长度示意图

根据能量原理已知M1=M2，因求得上述公式。再根据LEFM原理计算裂缝生成区段，公式如下：

(6)

式中，st—等效强度度值。

综合以上两组计算结果，得到有效裂缝扩展长度Δλ*：

(7)

在应变模型中，可参考单直线模型构建受力变形模型。横轴表示裂缝张开宽度，纵轴表示裂缝张开宽度为b时的闭合力；b0表示虚拟裂缝的最大张开度[10]。根据传统方法得出：虚拟裂缝最大张开度，与混凝土断裂能、抗拉强度以及最大骨料粒径有关。因此根据单直线模型，计算Δλ*的最终确定值：

Δλ*=Δλ-t0=4t0-t0=3t0

(8)

将矩形诱导缝的断裂问题，简化为无限大体深埋椭圆形裂缝问题，简化后的无限大体深埋椭圆裂缝模型如图2所示。

图2 无限大体深埋椭圆裂缝模型

图中的u、v、w，分别为模型的基本长度数据。确定有效裂缝扩展量后，再确认诱导缝等效强度，计算结果为：

(9)

式中，d—裂缝出现失稳扩展时的等效强度；Φ0—第二类完整椭圆积分；α—诱导缝四周，对裂缝应力强度因子的影响数值[11]。将公式两边同时除以混凝土大坝的抗拉强度d1，则得到

(10)

式中，λ=Δλ+Δλ*，计算结果为诱导缝的失稳相对等效强度。

同理计算另外两个阶段的等效强度，构建诱导缝等效强度双参数受力变形模型，以此建立双参数受力变形模型。

1.3 模拟混凝土大坝诱导缝变形过程

以常见的混凝土大坝弯曲切口梁试件为标准，模拟混凝土大坝中诱导缝变形过程。该模拟过程需要迭代计算，采用有限元法，先分析混凝土大坝的结构弹性，在裂纹尖端设置裂纹圆形弹性解析元，将普通有限元与其相连接，初步确定塑性区域尺寸数值。再将该值作为裂纹塑性初始值，分析混凝土大坝的结构塑性，该分析利用上一节构建的受力变形模型，将普通有限元与诱导缝圆形解析元相结合，根据裂纹尖端应力场与位移场，确定大坝在新的塑性区尺寸，并确定裂纹尖端张开或滑动位移量[12]。根据新的塑性区，重新计算数值，如此迭代计算结构塑性分析结果，直至裂纹尖端张开或滑动位移、塑性区尺寸完全收敛。混凝土大坝三点弯曲梁A1，在裂缝稳定扩展阶段的临界点，对应于最大轴向荷载时，裂缝尖端的切向应力计算结果，如图3(a)所示。

图3 大坝梁裂缝尖端切向应力分布图

根据图3(a)可知，临界裂缝稳定扩展阶段中，裂缝长度达到亚临界点时，裂缝就会发展到不稳定扩展阶段，此时的荷载虽然呈现下降趋势，但裂缝受先前作用力的影响，还是会继续扩展。而图3(b)是虚拟裂缝扩展长度达到F2时的裂缝尖端切向应力计算结果[13]。因此采用有限元计算方法，利用如图4所示的网格图，对混凝土大坝诱导缝变形过程进行数值模拟[14]。

图4 混凝土大坝弯曲梁网格剖分示意图

有限元法计算诱导缝受力变形数值，得到变形过程的全程模拟结果，如图5所示。

图5 大坝诱导缝变形过程全程模拟预期结果

根据图5曲线可以看出，需要有限元分析曲线，与实际值高度吻合，才能保证分析结果可靠，才能使构建的双参数受力变形模型，准确反馈裂纹稳定扩展和不稳定扩展的全过程[15]。至此利用上述方法，实现对混凝土大坝诱导缝受力变形的数值模拟。

2 实验研究

将此次提出的数值模拟方法作为实验A组，将传统方法[2- 3]，分别作为实验B组和实验C组，比较三个测试组的数值模拟效果。

选择1处大坝诱导缝受力变形区域记为N1，已知N1处存在5处变形节点，利用三个测试组，根据混凝土水力大坝基本数值，模拟大坝诱导缝受力变形过程。三个测试组的应力强度因子计算结果见表1—3。

表1 实验A组应力强度因子与解析解的比较

表2 实验B组应力强度因子与解析解的比较

表3 实验C组应力强度因子与解析解的比较

根据表1—3中，应力强度因子计算比较结果可知，实验A组的因子误差最小，而实验B组和实验C组的数值解高度相似，但其因子误差相对来说偏大，可见构建的模型会存在差异。对大坝诱导缝变形过程的全程模拟结果如图6所示。

图6 大坝诱导缝变形过程全程模拟的结果

根据图中的曲线可知，实验A组的有限元分析曲线，与混凝土水力大坝的实际曲线之间，没有极大的差异，可见模型得到的曲线与实际变化值高度吻合。实验B组和C组的曲线走势非常近似，但两个数值模拟方法构建的模型，受应力强度因子数值解的影响，生成的有限元分析曲线走势，与实际曲线存在微弱差异。利用下列方程，分别计算三个测试组数值模拟结果与实际曲线之间的相似度，结果见表4。

(11)

式中，e1—模型模拟诱导缝受力变形的有限元分析曲线数值；e2—混凝土大坝诱导缝受力变形的实际值；n—节点编号；i、j—节点坐标。

表4 模拟曲线与实际曲线之间的吻合程度

根据表中数据可知，实验A组的5个模拟节点数值，其相似度均在0.95以上，而实验B组和实验C组也有较高的相似度，但相对来说，此次研究的数值模拟方法效果更好。

3 结束语

社会经济高度发展下，建设的混凝土大坝规模会越来越大，能够服务的范围也会越来越广，因此其孔边裂缝、尖端裂缝还会出现，这些裂缝直接影响诱导缝受力变形程度，因此应力强度因子的计算是十分重要的。目前的应力强度因子计算方法，包括解析法、数值法以及实验法，这些方法各有特点，在今后的研究中，可以将上述计算方法进行融合，设计一个性能更加优越的计算方法，使应力强度因子的计算结果更加贴合真实值，进一步优化数值模拟效果。