CN111460568B - 一种混凝土重力坝运行期裂缝扩展判别方法 - Google Patents

Abstract

本发明公开了一种混凝土重力坝运行期裂缝扩展判别方法，包括步骤：确定大坝的几何、力学参数及不利工况；建立整体三维有限元模型，获取各种不利工况下大坝应力场和变形场；根据大坝裂缝所在位置和尺寸确定裂缝所在子区域；对每种不利荷载确定相应子区域的边界条件；重新优化裂缝子区域网格划分及节点布置；采用扩展有限元法建立非线性离散控制方程；求解非线性离散控制方程；判断裂缝是否扩展，若裂缝扩展则对扩展后的裂缝增加节点继续计算直至本工况裂缝不扩展，类似地再进行下一个工况直至所有不利工况全部模拟完。本发明实现了变结点扩展有限元和内聚力模型的深度融合，充分考虑了混凝土拟脆性材料的性质，裂缝扩展判别准确高效。

Description

一种混凝土重力坝运行期裂缝扩展判别方法

技术领域

本发明属于混凝土重力坝技术领域，具体涉及一种混凝土重力坝运行期裂缝扩展判别方法。

背景技术

我国现有大量的混凝土重力坝，其中许多重力坝存在深度裂缝。裂缝的存在严重影响大坝的整体性和安全性，尤其是沿上下游方向的深度裂缝。现有的裂缝分析方法将混凝土作为线弹性模型，没有准确考虑混凝土的拟脆性材料的性质。

发明内容

本发明的目的在于克服现有技术中的不足，提供一种混凝土重力坝运行期裂缝扩展判别方法，实现了扩展有限元、变结点单元和内聚力模型的深度融合，充分考虑了混凝土拟脆性材料的性质，裂缝扩展判别准确性高、效率高。

本发明提供了如下的技术方案：

一种混凝土重力坝运行期裂缝扩展判别方法，包括以下步骤：

收集混凝土重力坝的几何、力学参数，确定不利荷载组合；

建立大坝的整体三维有限元模型，对于每种不利荷载组合获取相应工况下整体应力场和变形场；

根据大坝裂缝所在位置和尺寸，明确裂缝所在区的子区域；

根据整体大坝的应力场和变形场确定相应不利工况下子区域的边界条件；

根据裂缝分布、子区域的几何材料和结构特征，以提高计算精度为目的重新优化细化裂缝子区域单元及网格划分；

采用扩展有限元法建立非线性离散控制方程；

采用最大拉应力准则模拟裂缝扩展；

求解非线性离散控制方程；

通过最大拉应力准则判断裂缝是否继续扩展，若扩展则根据裂缝扩展情况优化节点，再根据扩展有限元建立非线性离散控制方程模拟裂缝扩展直至本工况作用下裂缝不扩展为止，类似地再进行下一个工况直至所有不利工况全部模拟完。

进一步的，建立大坝包括基础的整体三维有限元模型，其中上下游方向及建基面以下深度各取两倍坝高，根据网格精度和计算稳定性要求进行网格剖分，通过不利荷载组合采用混凝土徐变本构模型和多场耦合有限单元法获得整体大坝的应力场和变形场。

进一步的，根据大坝裂缝所在位置和尺寸，截取裂缝所在平面的包括裂缝长度的多边形区域作为覆盖裂缝的子空间：如果裂缝位于断面边缘，则向坝体内部深度取3倍最大缝长，高度上下各取1.5倍最大缝长形成一个矩形子区域；如果裂缝位于断面中部，又分两种情况：对于连通裂缝区域，裂缝的端部边界到子区域相应边界直线距离不小于整个缝长，且子区域长宽比不大于2；对于邻近的分区裂缝，采用多个连在一起的满足上述要求的多个连在一起的矩形进行完整覆盖。

进一步的，确定不利工况下子区域的边界条件是根据整体大坝的应力场和变形场采用子结构法确定裂缝子区域的边界条件。

进一步的，优化细化裂缝子区域单元及网格划分的方法包括以下步骤：

对于四结点四边形单元，初始裂尖单元的非连续位移场表示如下：

式中，N^k是裂尖单元上的结点集，b_i是结点自由度，N_i(x^*)是四结点四边形单元内点x^*的标准结点形函数值，H(x)是修正的Heaviside阶跃函数，在裂缝上方取+1，下方取-1；

开裂单元的位移场是连续与非连续位移场的和，借助四结点四边形单元将内聚裂缝模拟的位移逼近写为：

式中，N^s是离散域中所有结点的集合，N^cut是被裂缝完全切割单元的结点集，u_i和a_i分别是结点自由度向量和结点加强变量，N_i(x)是结点形函数；

由式(15)可知，裂缝张开位移可表示为：

采用变结点单元连接两个不同尺度的网格区域，基于具有特殊基底的一般点插值的概念而建立变结点单元的形函数，这些基底在二维域中具有不共线性，在此框架下，位移逼近u^h(x)表示为：

式中，N_p是插值中的样点个数，N_i是第i个点的形函数矩阵，u_i是变结点向量，a^T是未知系数矩阵，p(ξ)是多项式基底的列向量；

对于(4+k+m)结点单元，其中k是单元顶部和底部边缘的额外结点数，m是左侧和右侧边缘的额外结点数，多项式基底如下：

式中η和ξ描述了等参单元中的局部坐标；

点插值如下：

u^h(ξ)＝a^Tp(ξ)＝U^Tq^-1p(ξ) (19)

其中(4+k+m)结点单元的形状函数定义为：

[N₁,…,N_4+k+m]^T＝q^-1p(ξ) (20)

其中

q＝[p(ξ₁),…,p(ξ_4+k+m)] (21)

变结点单元保留了任意两个相邻结点之间的线性插值，因此它们可以自然地连接不同尺度网格的单元。

进一步的，建立非线性离散控制方程时得到的方程组为：

Kδ＝R (25)

式中，K和R分别是整体刚度矩阵与外结点力向量，R由大坝整体有限元子空间方法获得，δ是结点未知量向量；

每个单元的局部刚度矩阵给出如下：

且

其中形函数导数定义为：

每个单元的外结点力表示为：

其中：

其中，

和

分别为常规单元荷载列向量和被裂纹穿过单元荷载列向量。

进一步的，模拟裂缝扩展的步骤包括：

当虚拟裂缝尖端应力张量在法向上的投影达到材料的抗拉强度时，裂缝将发生扩展，如下所示：

n₂·C·B(x_t)δ＝S·δ＝f_t (33)

其中

n₂＝[l² m² 2lm] (34)

S＝n₂·C·B(x_t) (35)

其中l和m是局部坐标系中虚拟裂缝尖端法向量n的分量，C为弹性矩阵，x_t是虚拟裂缝尖端的位置。

此外，裂缝扩展方向是由最大环向应力准则来确定的：

其中，

和

分别是在外力作用下，模式I和模式II的应力强度因子。

进一步的，求解非线性离散控制方程的步骤包括：

采用弧长法迭代求解式(25)的离散控制方程。荷载因子λ和结点变量δ视为未知量，不平衡力

计算公式为：

其中，

和f^coh(δ)分别为外力τ₀和内聚力τ^c的等效结点荷载向量；

在第i次迭代时，使用泰勒展开式对不平衡力

进行一阶展开：

令

得到：

将式(39)代入δⁱ⁺¹＝δⁱ+Δδⁱ，并与式(33)联立，得到：

荷载因子增量此处定义为：

未知增量是由方程(39)～(41)推导出的，然后第(i+1)次迭代的荷载和位移由λⁱ⁺¹＝λⁱ+Δλⁱ和δⁱ⁺¹＝δⁱ+Δδⁱ给出；

当目前的不平衡力范数与前一步范数之比小于一定值时，认为满足平衡条件，即解是收敛的。

与现有技术相比，本发明的有益效果是：本发明通过收集混凝土重力坝的几何、力学参数，大坝整体数值建模，明确裂缝区子区域，确定相应不利工况下子区域的边界条件，优化裂缝子区域单元划分，建立非线性离散控制方程，模拟裂缝扩展，求解非线性控制方程，最终完成裂缝扩展判断，实现了扩展有限元、变结点单元和内聚力模型的深度融合，充分考虑了混凝土拟脆性材料的性质，裂缝扩展判别准确性高、效率高。

附图说明

图1是混凝土重力坝运行期裂缝扩展判别方法的流程示意图；

图2是实施例中二维平面开裂体的结构示意图；

图3是实施例中所采用的线性内聚定律的示意图；

图4是连通裂缝子区域示意图；

图5是分区裂缝子区域示意图；

图6是实施例中四结点四边形单元的示意图；

图7是图6等参单元1-2-3-4中结点5、6的局部坐标示意图；

图8是实施例中两种不同尺度单元的网格示意图；

图9是实施例中(4+k+m)结点单元示意图。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

如图2所示的二维平面开裂体Ω以及边界Γ＝Γ_u+Γ_t+Γ_c，

其中Γ_u、Γ_t和Γ_c分别是位移边界、应力边界和裂缝面，不考虑体积力，边界条件描述如下：

σ·n＝λτ₀,onΓ_t (2)

σ⁺·n⁺＝σ^-·n^-＝τ^c,onΓ_coh (4)

式中

是梯度算子，σ是柯西应力张量，u是位移，n是边界Γ上的向外单位法向量，

是规定位移，τ₀是标准牵引力，λ是荷载因子，并且n⁺＝-n^-＝n^c是裂缝面上的向外单位法向量，σ^±和τ^c是裂缝面上的应力和法向牵引力，Γ_coh是内聚区裂缝面。

考虑小变形，以下关系成立：

式中ε是应变张量；

加之本构方程：

σ＝Cε,inΩ (6)

τ^c＝τ^c(w),onΓ_coh (7)

其中C是四阶弹性张量的矩阵形式，w是裂缝开口，定义如下：

w＝u^--u⁺,onΓ_coh (8)

u^-和u⁺分别是裂缝表面的位移。

将虚功原理应用于等式(6)，平衡方程的弱形式可以获得如下：

其中v是容许虚位移，v^±是裂缝面上的容许虚位移。

考虑式(4)中的连续性条件和w(v)＝v^--v⁺，式(9)可以被改写为：

一般来说，裂缝面上的内聚力τ^c是法向和切向裂缝张开量的函数，为了简化和充分应用混合内聚力模型，本发明只考虑法向牵引模式，而忽略了剪切牵引，这样更有利于对裂缝的准确模拟。本发明采用线性内聚力进行分析，如图3，法向内聚力τ^c表示为：

其中f_t是抗拉强度，w和w₀分别是法向裂缝张开位移和临界裂缝张开位移，内聚断裂能量定义为

如图1所示，一种混凝土重力坝运行期裂缝扩展判别方法，包括以下步骤：

步骤一、收集混凝土重力坝的几何、力学参数，确定不利荷载组合

收集重力坝和基础一定范围内的结构和几何尺寸、材料力学参数、材料热力学参数以及极端情况下的水荷载和温度荷载，通过Couple方法确定不利荷载组合。

步骤二、大坝整体数值建模

建立大坝包括基础的整体三维有限元模型，其中上下游方向及建基面以下深度各取两倍坝高，根据网格精度和计算稳定性要求进行网格剖分，通过不利荷载组合采用混凝土徐变本构模型和多场耦合有限单元法获得整体大坝的应力场和变形场。

步骤三、明确裂缝所在区的子区域

根据大坝裂缝所在位置和尺寸，截取裂缝所在平面的包括裂缝长度的多边形区域作为覆盖裂缝的子空间：如果裂缝位于断面边缘，则向坝体内部深度取3倍最大缝长，高度上下各取1.5倍最大缝长形成一个矩形子区域；如果裂缝位于断面中部，又分两种情况：对于连通裂缝区域，裂缝的端部边界到子区域相应边界直线距离不小于整个缝长，且子区域长宽比不大于2，如图4为连通裂缝子区域示意图；对于邻近的分区裂缝，采用多个连在一起的满足上述要求的多个连在一起的矩形进行完整覆盖，如图5为分区裂缝子区域示意图。

步骤四、确定相应不利工况下子区域的边界条件

根据步骤二计算结构采用子结构法确定步骤三裂缝子空间的边界条件。

步骤五、重新优化细化裂缝子区域单元及网格划分

根据下列方法对裂缝子空间进行单元重新划分：

如图4的四结点四边形单元，1-4为裂尖切割边，要在边1-4上构建富集场，过p点做辅助线交单元边界于5、6两点，裂尖富集单元即为单元1-5-6-4，结点1、4被富集了，初始裂尖单元的非连续位移场表示如下：

u_disc＝b₁N₁(x^*)(H(x^*)-H(x₁))+b₄N₄(x^*)(H(x^*)-H(x₄)) (12)

其中N_i(x^*)是单元1-5-6-4内的点x^*的标准结点形函数值，b_i是结点自由度，H(x)是修正的Heaviside阶跃函数，在裂缝上方取+1，下方取-1；如图4所示，等参单元1-2-3-4中结点5、6的局部坐标如图5所示。

相似的，如果裂缝切割边界为1-2，初始裂尖单元的非连续位移场则表示为：

u_disc＝b₁N₁(x^*)(H(x^*)-H(x₁))+b₂N₂(x^*)(H(x^*)-H(x₂)) (13)

则初始裂尖单元的非连续位移场表示如下：

式中，N^k是裂尖单元上的结点集，N_i(x^*)是四结点四边形单元内点x^*的标准结点形函数值；

开裂单元的位移场是连续与非连续位移场的和，借助四结点四边形单元将内聚裂缝模拟的位移逼近写为：

式中，N^s是离散域中所有结点的集合，N^cut是被裂缝完全切割单元的结点集，u_i和a_i分别是结点自由度向量和结点加强变量，N_i(x)是结点形函数；

由式(15)可知，裂缝张开位移可表示为：

本实施例采用多尺度网格提高计算准确性，即只在关键区域(预期开裂)进行精细网格化，而粗尺度网格化运用于其余部分，在两种不同尺度单元之间存在一层变结点单元，如图6中的阴影部分所示；

本实施例采用变结点单元连接两个不同尺度的网格区域，基于具有特殊基底的一般点插值的概念而建立变结点单元的形函数，这些基底在二维域中具有不共线性，在此框架下，位移逼近u^h(x)表示为：

式中，N_p是插值中的样点个数，N_i是第i个点的形函数矩阵，u_i是变结点向量，a^T是未知系数矩阵，p(ξ)是多项式基底的列向量；

如图7显示了一个(4+k+m)结点单元，其中k是单元顶部和底部边缘的额外结点数，m是左侧和右侧边缘的额外结点数，多项式基底如下：

式中η和ξ描述了等参单元中的局部坐标；

点插值如下：

u^h(ξ)＝a^Tp(ξ)＝U^Tq^-1p(ξ) (19)

其中(4+k+m)结点单元的形状函数定义为：

[N₁,…,N_4+k+m]^T＝q^-1p(ξ) (20)

其中

q＝[p(ξ₁),…,p(ξ_4+k+m)] (21)

变结点单元保留了任意两个相邻结点之间的线性插值，因此它们可以自然地连接不同尺度网格的单元。

步骤六、采用扩展有限元法建立非线性离散控制方程

将式(10)离散得到以下方程组：

Kδ＝R (25)

式中，K和R分别是整体刚度矩阵与外结点力向量，R由大坝整体有限元子空间方法获得，δ是结点未知量向量；

每个单元的局部刚度矩阵给出如下：

且

其中形函数导数定义为：

每个单元的外结点力表示为：

其中：

其中，

和

分别为常规单元荷载列向量和被裂纹穿过单元荷载列向量。

步骤七、模拟裂缝扩展

模拟裂缝扩展采用最大拉应力准则，即当虚拟裂缝尖端应力张量在法向上的投影达到材料的抗拉强度时，裂缝将发生扩展，如下所示：

n₂·C·B(x_t)δ＝S·δ＝f_t (33)

其中

n₂＝[l² m² 2lm] (34)

S＝n₂·C·B(x_t) (35)

其中l和m是局部坐标系中虚拟裂缝尖端法向量n的分量，C为弹性矩阵，x_t是虚拟裂缝尖端的位置。

此外，裂缝扩展方向是由最大环向应力准则来确定的：

其中，

和

分别是在外力作用下，模式I和模式II的应力强度因子。

步骤八、求解非线性控制方程

采用弧长法迭代求解式(25)的离散控制方程。荷载因子λ和结点变量δ视为未知量，不平衡力

计算公式为：

其中，

和f^coh(δ)分别为外力τ₀和内聚力τ^c的等效结点荷载向量；

在第i次迭代时，使用泰勒展开式对不平衡力

进行一阶展开：

令

得到：

将式(39)代入δⁱ⁺¹＝δⁱ+Δδⁱ，并与式(33)联立，得到：

荷载因子增量此处定义为：

未知增量是由方程(39)～(41)推导出的，然后第(i+1)次迭代的荷载和位移由λⁱ⁺¹＝λⁱ+Δλⁱ和δⁱ⁺¹＝δⁱ+Δδⁱ给出；

当目前的不平衡力范数与前一步范数之比小于一定值时，认为满足平衡条件，即解是收敛的。

步骤九、裂缝扩展判断

通过最大拉应力准则判断裂缝是否继续扩展，若扩展则根据裂缝扩展情况优化节点，再根据扩展有限元建立非线性离散控制方程模拟裂缝扩展直至本工况作用下裂缝不扩展为止，类似地再进行下一个工况直至所有不利工况全部模拟完。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims (8)

1.一种混凝土重力坝运行期裂缝扩展判别方法，其特征在于，包括以下步骤：

收集混凝土重力坝的几何、力学参数，确定不利荷载组合；

建立大坝的整体三维有限元模型，对于每种不利荷载组合获取相应工况下整体应力场和变形场；

根据大坝裂缝所在位置和尺寸，明确裂缝所在区的子区域；

根据整体大坝的应力场和变形场确定相应不利工况下子区域的边界条件；

根据裂缝分布、子区域的几何材料和结构特征，以提高计算精度为目的重新优化细化裂缝子区域单元及网格划分；

采用扩展有限元法建立非线性离散控制方程；

采用最大拉应力准则模拟裂缝扩展；

求解非线性离散控制方程；

通过最大拉应力准则判断裂缝是否继续扩展，若扩展则根据裂缝扩展情况优化节点，再根据扩展有限元建立非线性离散控制方程模拟裂缝扩展直至本工况作用下裂缝不扩展为止，类似地再进行下一个工况直至所有不利工况全部模拟完。

2.根据权利要求1所述的混凝土重力坝运行期裂缝扩展判别方法，其特征在于，建立大坝包括基础的整体三维有限元模型，其中上下游方向及建基面以下深度各取两倍坝高，根据网格精度和计算稳定性要求进行网格剖分，通过不利荷载组合采用混凝土徐变本构模型和多场耦合有限单元法获得整体大坝的应力场和变形场。

3.根据权利要求1所述的混凝土重力坝运行期裂缝扩展判别方法，其特征在于，根据大坝裂缝所在位置和尺寸，截取裂缝所在平面的包括裂缝长度的多边形区域作为覆盖裂缝的子空间：如果裂缝位于断面边缘，则向坝体内部深度取3倍最大缝长，高度上下各取1.5倍最大缝长形成一个矩形子区域；如果裂缝位于断面中部，又分两种情况：对于连通裂缝区域，裂缝的端部边界到子区域相应边界直线距离不小于整个缝长，且子区域长宽比不大于2；对于邻近的分区裂缝，采用多个连在一起的满足上述要求的矩形进行完整覆盖。

4.根据权利要求1所述的混凝土重力坝运行期裂缝扩展判别方法，其特征在于，确定不利工况下子区域的边界条件是根据整体大坝的应力场和变形场采用子结构法确定裂缝子区域的边界条件。

5.根据权利要求1所述的混凝土重力坝运行期裂缝扩展判别方法，其特征在于，优化细化裂缝子区域单元及网格划分的方法包括以下步骤：

对于四结点四边形单元，初始裂尖单元的非连续位移场表示如下：

式中，N^k是裂尖单元上的结点集，b_i是结点自由度，N_i(x^*)是四结点四边形单元内点x^*的标准结点形函数值，H(x)是修正的Heaviside阶跃函数，在裂缝上方取+1，下方取-1；

开裂单元的位移场是连续与非连续位移场的和，借助四结点四边形单元将内聚裂缝模拟的位移逼近写为：

式中，N^s是离散域中所有结点的集合，N^cut是被裂缝完全切割单元的结点集，u_i和a_i分别是结点自由度向量和结点加强变量，N_i(x)是结点形函数；

由式(15)可知，裂缝张开位移可表示为：

采用变结点单元连接两个不同尺度的网格区域，基于具有特殊基底的一般点插值的概念而建立变结点单元的形函数，这些基底在二维域中具有不共线性，在此框架下，位移逼近u^h(x)表示为：

式中，N_p是插值中的样点个数，N_i是第i个点的形函数矩阵，u_i是变结点向量，a^T是未知系数矩阵，p(ξ)是多项式基底的列向量；

对于(4+k+m)结点单元，其中k是单元顶部和底部边缘的额外结点数，m是左侧和右侧边缘的额外结点数，多项式基底如下：

式中η和ξ描述了等参单元中的局部坐标；

点插值如下：

u^h(ξ)＝a^Tp(ξ)＝U^Tq^-1p(ξ) (19)

其中(4+k+m)结点单元的形状函数定义为：

[N₁,...,N_4+k+m]^T＝q^-1p(ξ) (20)

其中

q＝[p(ξ₁),...,p(ξ_4+k+m)] (21)

变结点单元保留了任意两个相邻结点之间的线性插值，因此它们可以自然地连接不同尺度网格的单元。

6.根据权利要求1所述的混凝土重力坝运行期裂缝扩展判别方法，其特征在于，建立非线性离散控制方程时得到的方程组为：

Kδ＝R (25)

式中，K和R分别是整体刚度矩阵与外结点力向量，R由大坝整体有限元子空间方法获得，δ是结点未知量向量；

每个单元的局部刚度矩阵给出如下：

且

其中形函数导数定义为：

每个单元的外结点力表示为：

r_i＝[r_i ^u r_i ^a]^T (30)

其中

其中，r_i ^u和r_i ^a分别为常规单元荷载列向量和被裂纹穿过单元荷载列向量。

7.根据权利要求1所述的混凝土重力坝运行期裂缝扩展判别方法，其特征在于，模拟裂缝扩展的步骤包括：

当虚拟裂缝尖端应力张量在法向上的投影达到材料的抗拉强度时，裂缝将发生扩展，如下所示：

n₂·C·B(x_t)δ＝S·δ＝f_t (33)

其中

n₂＝[l² m² 2lm] (34)

S＝n₂·C·B(x_t) (35)

其中l和m是局部坐标系中虚拟裂缝尖端法向量n的分量，C为弹性矩阵，x_t是虚拟裂缝尖端的位置。

此外，裂缝扩展方向是由最大环向应力准则来确定的：

其中，

和

分别是在外力作用下，模式I和模式II的应力强度因子。

8.根据权利要求1所述的混凝土重力坝运行期裂缝扩展判别方法，其特征在于，求解非线性离散控制方程的步骤包括：

采用弧长法迭代求解式(25)的离散控制方程。荷载因子λ和结点变量δ视为未知量，不平衡力

计算公式为：

其中，

和f^coh(δ)分别为外力τ₀和内聚力τ^c的等效结点荷载向量；

在第i次迭代时，使用泰勒展开式对不平衡力

进行一阶展开：

令

得到：

将式(39)代入δⁱ⁺¹＝δⁱ+Δδⁱ，并与式(33)联立，得到：

荷载因子增量此处定义为：

未知增量是由方程(39)～(41)推导出的，然后第(i+1)次迭代的荷载和位移由λⁱ⁺¹＝λⁱ+Δλⁱ和δⁱ⁺¹＝δⁱ+Δδⁱ给出；

当目前的不平衡力范数与前一步范数之比小于一定值时，认为满足平衡条件，即解是收敛的。

Images