CN110765695B - 一种基于高阶有限元法获取混凝土重力坝裂纹扩展路径的模拟计算方法 - Google Patents

Abstract

本发明涉及一种基于高阶有限元法获取混凝土重力坝裂纹扩展路径的模拟计算方法，属于水利水电工程领域，特别是对已含裂纹的混凝土重力坝裂纹扩展的路径进行分析计算。利用高阶有限元法计算出具有较高精度的混凝土重力坝应力场和位移场，进而得出高精度的应力强度因子及相当应力强度因子，最后获得相对准确的混凝土重力坝裂纹扩展路径。本计算方法能够直接应用于混凝土重力坝的工程实际问题中，对含初始微裂纹的混凝土重力坝进行裂纹扩展路径的模拟计算，并且计算结果较现有其他计算方法精度相对较高。

Description

一种基于高阶有限元法获取混凝土重力坝裂纹扩展路径的模 拟计算方法

技术领域

本发明涉及一种基于高阶有限元法获取混凝土重力坝裂纹扩展路径的模拟计算方法，属于水利水电工程领域，特别是对已含裂纹的混凝土重力坝裂纹扩展的路径进行分析计算。

背景技术

近些年来，我国为了防洪发电以及灌溉修建了许多大大小小的水电站及水库，而作为对水电站与水库极为重要的挡水建筑也受到了极大的关注。混凝土重力坝由于其适用性、可靠性、稳定性及施工的便捷性在挡水建筑物中占用举足轻重的地位。由于混凝土材料自身的抗拉性能较弱以及混凝土重力坝的大体积特性，以致混凝土重力坝会不可避免地在施工中产生微裂缝。初始微裂缝使混凝土重力坝的安全受到潜在威胁，并且随着混凝土重力坝运行时间的增加，初始微裂缝可能会进一步扩展，最终可能会造成非常严重的事故。因此在裂纹问题方面，混凝土重力坝安全评价的首要问题是判断已有裂纹在什么情况下会扩展，以及预判在各种可能工况下的裂纹扩展路径及其对应的安全程度，以便在裂纹进一步扩展危及坝体前采取控制防护措施。但是，目前对裂纹扩展方面的研究并不多。

现在关于不连续问题多采用扩展有限元方法进行分析，在裂纹扩展方面扩展有限元方法也有其独特的优势。但扩展有限元方法是基于传统有限元方法发展起来的，它具有精度较差，收敛速度较慢的缺点。本发明涉及的基于高阶有限元法获取混凝土重力坝裂纹扩展路径的方法利用了高阶有限元法高精度，收敛速度快以及稳定性较好的优点，从而获得相对准确的混凝土重力坝裂纹扩展路径。

使用扩展有限元方法对裂纹扩展问题进行研究的几乎都是基于传统有限元方法，因此存在如下缺点：

(1)传统有限元方法对裂纹扩展问题缺少有效的误差估计模式，最终计算结果精度的控制比较依赖于研究者的经验。为了提高计算结果的精度，需要进行多次重新划分网格，以提高有限元解的收敛性；

(2)传统有限元方法基于低阶插值函数，因此在模拟裂纹尖端附近的高梯度应力和应变场时存在不足，为了获得较高的计算精度，需要进行多次加密网格以确保足够的网格密度，前处理与计算成本较高；

(3)由于裂尖附近应力场的奇异性，基于传统有限元方法计算得到的裂尖附近应力与应变场的误差较大，进而由其导出的应力强度因子精度较低；若想要获得较高精度的应力强度因子，则前处理和计算成本将随之提高。

专利申请号为201910063707X的发明专利主要利用了p型有限元法结合围线积分法来导出复合型应力强度因子，而本发明选择公开基于高阶有限元法通过循环计算得到满足工程精度要求的混凝土重力坝裂纹扩展路径。

本发明基于国家自然科学基金资助项目(资助编号：51769011)，提出了一种将高阶有限元法应用于水利水电工程领域，从而获取混凝土重力坝裂纹扩展路径的新型计算技术。

发明内容

针对现有混凝土重力坝裂纹扩展路径模拟计算方法的缺点及不足，本发明提供一种基于高阶有限元法获取混凝土重力坝裂纹扩展路径的模拟计算方法。本发明通过以下技术方案实现。

一种基于高阶有限元法获取混凝土重力坝裂纹扩展路径的模拟计算方法，利用高阶有限元法计算出具有较高精度的混凝土重力坝应力场和位移场，进而得出高精度的应力强度因子及相当应力强度因子，最后获得相对准确的混凝土重力坝裂纹扩展路径，包括以下步骤：

步骤1、确定混凝土重力坝工程模型的几何尺寸，并创建混凝土重力坝的有限元模型；

步骤2、确定初始裂纹的位置及初始裂纹的长度，并将初始裂纹创建在混凝土重力坝有限元模型上；

步骤3、设置初始开裂方向、初始开裂位置及开裂步长；

步骤4、设置插值多项式阶数及循环次数；

每一次循环自动划分网格来进行混凝土重力坝裂纹扩展路径分析，得到每一循环步裂尖的应力强度因子，并求出相当应力强度因子，将相当应力强度因子和混凝土断裂韧度相比较，相当应力强度因子大于断裂韧度的循环步为有效循环，其余循环步舍去，判断有效循环步最后一步相当应力强度因子与混凝土断裂韧度之差是否在误差允许范围内，若不在误差允许范围内，则返回步骤3调小开裂步长并增大循环次数，继续进行计算。

所述步骤1包括混凝土重力坝几何模型的建立，混凝土材料参数的设置，载荷及边界条件的施加。

所述步骤4具体包括以下几步：

步骤4.1、对插值多项式阶数和循环次数进行假定，其中插值多项式阶数为1～8；

步骤4.2、通过确定了插值多项式阶数的高阶有限元法计算得到含裂纹混凝土重力坝第一次循环的位移场和应力应变场，并且得到裂纹扩展角；

步骤4.3、利用围线积分法导出每一次循环的应力强度因子K_I和K_II；

步骤4.4、计算得到第一次循环的相当应力强度因子

步骤4.5、通过裂纹扩展角和裂纹扩展步长完成裂纹扩展第一步循环；

步骤4.6、开始第二步循环，重复步骤4.2至步骤4.6，完成裂纹扩展第二步循环，以此类推，直到完成所有循环次数，将所有计算得到的每一步循环的相当应力强度因子依次与混凝土断裂韧度对比；

若所有循环步的相当应力强度因子均大于混凝土断裂韧度对比，则增大循环次数，返回步骤4.1，若存在某些循环步的相当应力强度因子小于混凝土断裂韧度，则找出最接近且大于混凝土断裂韧度的那一步的相当应力强度因子，求得该步相当应力强度因子和混凝土断裂韧度的相对误差，看相对误差是否满足工程精度要求，如果满足，则裂纹扩展完成，如果不满足且相对误差与工程精度要求相差较小，则可尝试增大插值多项式阶数重复上述计算，如果不满足且相对误差与工程精度要求相差较大，则减小裂纹扩展步长并增大循环次数，减小后的裂纹扩展步长与增大后的循环次数的乘积应大于相当应力强度因子最接近且大于混凝土断裂韧度的那一步的循环步数与减小前的裂纹扩展步长的乘积，再重复上述步骤，直至满足工程精度要求。

上述步骤4.2中计算每一次循环的位移场和应力应变场具体步骤如下：

步骤4.2.1、依据含初始裂纹的混凝土重力坝有限元模型自动划分网格，结合所指定的插值多项式阶数，求解结构位移列阵a

依据方程：Ka＝F (1)

其中：K为总体整体刚度矩阵，K＝∑_eG^TK_eG，(K_e为单元刚度矩阵)；

a为结构位移列阵；

F为结构节点载荷列阵，F＝∑_eG^Tf_e，(e表示单元的)；

G表示结构节点自由度与单元节点自由度的转换矩阵，T表示转置；

再依据单元刚度矩阵的公式K_e＝∫_ΩB^TDBdΩ (2)

且单元等效节点载荷列阵

其中：单元内部节点力

外部节点力

上式(2)至上式(5)中，Ω表示为单元内部；B＝LN_I，B表示应变矩阵，L表示微分算子，N_I表示形函数矩阵，I表示单元的标号，Γ_t表示单元的外部边界；D表示应力矩阵，b表示体力，

表示荷载边界条件，这些量由步骤1中确定的混凝土材料属性、载荷及边界条件计算得到；

把上式(2)至上式(5)代入到上式(1)中，再对线性方程组(1)进行求解，则得到了结构位移列阵a；再依据公式a_I＝Ga进一步得到了单元位移列阵a_I；最后依据公式u＝N_Ia_I，ε＝Lu，σ＝Dε＝DBa_I，求得位移场u，应变场ε及应力场σ。

其中形函数矩阵N_I的计算采用高阶有限元法中的高阶插值多项式，典型的二维高阶有限元法形函数以勒让德正交多项式为基，以二维的四边形单元(如图1所示)为例说明高阶有限元法，在坐标系ξOη下，点P₁(-1,-1)、P₂(-1,1)、P₃(1,1)和P₄(-1,1)分别为四边形单元的四个顶点，Γ₁(η＝-1)、Γ₂(ξ＝1)、Γ₃(η＝1)和Γ₄(ξ＝-1)分别为四边形单元的四条边，构造形函数的步骤如下：

①当p＝1时，节点模式形状函数与常规的拉格朗日型插值函数相同，以四节点四边形单元为例展开如下：

②当p≥2时，边模式形状函数N_i ^(k)是与对应边Γ_k相关联的形状函数，与边

Γ₁＝{(ξ,-1),-1≤ξ≤1}相关联的形状函数定义如下：

其中：

P_n(t)为阶数n≥0的勒让德多项式

类似地，定义与边Γ_k(2≤k≤4)相关联的形状函数如下：

③当p≥4时，内部模式形状函数：

本发明采用高阶有限元法，其插值多项式的阶数可以从p＝1逐渐增大，当提升插值多项式的阶数时，先前低阶的刚度矩阵可以继续使用，只需要增加高阶部分的计算，避免了重复计算低阶部分的刚度矩阵，具有较为良好的沿用性，减少了前处理的计算，节约了计算成本。

步骤4.2.2、依据所求得的应力场和位移场，通过围线积分法导出应力强度因子，进而求出裂纹扩展角并结合确定的裂纹扩展步长完成这一步裂纹扩展；同时求出相当应力强度因子用以判断该步裂纹扩展是否为有效步；

本发明利用了围线积分法来导出应力强度因子。围线积分法是基于Betti功互等定理的一种方法，采用已求得的积分路径上的位移与应力应变导出混凝土重力坝裂纹的应力强度因子，并且该方法在导出应力强度因子时是超收敛的。如图2所示的二维弹性裂纹场Ω，在裂纹前端定义了局部笛卡尔坐标系(x₁,x₂)，记u(x₁,x₂)为局部坐标系(x₁,x₂)下弹性裂纹场Ω中的位移场。

表示在局部极坐标系(r,θ)下包围裂纹前端的邻域Ω_s中的位移场。当子域Ω_s已经足够接近裂纹尖端时，假设导出区域完全处在平面应变状态，则子域Ω_s内部的位移场可做如下近似：

式中，K_I与K_II分别对应Ⅰ型与Ⅱ型应力强度因子，G是剪切模量，κ是Kolosov常数,

以下用T^(u)(r,θ)来表示从位移场中导出的牵引力向量。应力强度因子K_I与K_II可以通过围线Γ₂，Γ₃及Γ₄上的线积分计算获得：

式中：

分别表示Ⅰ型与Ⅱ型应力强度因子的导出函数。

与

分别表示由

与

导出的牵引力向量。

式(6)与(7)算得的结果是精确的，但u的解及相应牵引力T^(u)是未知的，本发明采用高阶有限元法获得一个近似的数值解。围线Γ₂，Γ₃及Γ₄独立于用高阶有限元法求解近似数值解的网格，并且围线Γ₂的选取是较为灵活的，不必靠近裂纹的前端。由于裂纹前端的最内层单元中数值解结果的精度较低，此处求得的应力强度因子的误差较大，因此在选择围线Γ₂时，一般选择包围住裂尖附近最内层的单元。所述步骤4.2.2中围线积分法的围线远离并且包含裂尖，而且由于裂尖的外层单元受到奇异性的影响，所以围线也应当包含裂尖外的最内层单元。

上述4.2.1中混凝土重力坝有限元模型是自动划分网格，而无需提前划分网格或对网格尺寸进行设置，每一步循环自动划分网格主要是以前一步完成扩展的裂纹为单元边界，从而使对混凝土重力坝裂纹性态的描述更加简便。

本发明的有益效果是：

(1)本方法能够降低计算成本，加快收敛速度，提高计算精度及计算效率。

(2)本计算方法能够直接应用于混凝土重力坝的工程实际问题中，对含初始微裂纹的混凝土重力坝进行裂纹扩展路径的模拟计算，并且计算结果较现有其他计算方法精度相对较高。

附图说明

图1是本发明涉及的标准四边形母单元示意图；

图2是本发明涉及的全局坐标及局部坐标系下裂纹尖端邻域的示意图；

图3是本发明实施例1中含初始裂纹的混凝土重力坝的几何模型图；

图4是本发明实施例1中混凝土重力坝有限元模型计算得到的位移场云图；

图5是本发明实施例1中混凝土重力坝有限元模型计算得到的应力场云图；

图6是本发明实施例1中混凝土重力坝裂纹扩展路径示意图；

图7是本发明实施例1中混凝土重力坝变形图；

图8是本发明实施例1中扩展完成后的混凝土重力坝裂尖局部应力云图。

具体实施方式

下面结合附图和具体实施方式，对本发明作进一步说明。

实施例1

含初始裂纹的混凝土重力坝几何模型，如3所示其中坝高H＝100m，坝顶宽A＝6m，坝底宽B＝20，上游坝坡垂直，下游坝坡坡度为1:0.7；坝体为混凝土材料，考虑混凝土材料弹性模量E＝22GPa，泊松比ν＝0.250，密度ρ＝2700kg/m³；混凝土的断裂韧度K_c＝575314Pa·m¹/₂。坝体上游面承受静水压力，考虑最不利荷载情况，假定水位与坝顶高程齐平。假定在坝体上游面下部位置存在一条裂纹，且裂纹长度L＝10.00m，为方便描述与计算，考虑以上游面与坝底面交线的坝踵处作为坐标原点建立坐标系xoy，则裂纹两端端点的坐标为(0.00,13.5)和(10.00,13.5)。

步骤1、确定混凝土重力坝工程模型的几何尺寸，并创建混凝土重力坝的有限元模型。根据混凝土重力坝的几何尺寸建立几何模型；根据确定的材料属性参数建立有限元模型的材料类型，将相应的材料类型(E＝22GPa,ν＝0.250)赋予相对应的区域；根据确定的载荷和边界条件，在相应的边界上施加相对应的载荷(与坝顶齐平上游水位的上游水压力)和边界条件(坝底固定位移约束)，由于坝体自重对应力强度因子影响较小，因此本例忽略坝体自身重力；

步骤2、确定初始裂纹的位置(0,13.5)及初始裂纹的长度L＝10m，并将初始裂纹((0,13.5)-(10,13.5))创建在混凝土重力坝有限元模型上；

步骤3、设置初始开裂方向0°、初始开裂位置(10,13.5)及开裂步长a＝0(a为裂纹向左边开裂的开裂步长，因为此裂纹为边缘裂纹，因此为零),b＝0.05(b为裂纹向右边开裂的开裂步长，经多次重复计算，为满足精度要求故取b＝0.05)；

步骤4、经过多次试算，为满足工程精度要求，选择设置插值多项式阶数p＝5及循环次数为725次；

计算每一次循环的位移场和应力应变场，插值多项式阶数p＝5、第715次循环得到的位移场云图如图4所示，通过计算所得的应力场，应力场云图如图6所示

每一次循环自动划分网格来进行混凝土重力坝裂纹扩展路径分析，得到每一循环步裂尖的应力强度因子，并求出相当应力强度因子(考虑相当应力强度因子为

将相当应力强度因子和混凝土断裂韧度相比较，相当应力强度因子大于断裂韧度的循环步为有效循环，其余循环步舍去，判断有效循环步最后一步相当应力强度因子与混凝土断裂韧度之差是否在误差允许范围内，若不在误差允许范围内，则返回步骤3、4调小开裂步长并增大循环次数，继续进行计算。

经上述步骤，计算得裂纹扩展路径上裂尖点的坐标表如下(由于篇幅问题，本例只给出了725次循环中较为重要的702次-715次循环，716次-725次循环为无效循环)：

表1裂纹扩展路径上裂尖点的坐标及判定裂纹是否扩展表

本例第716次循环步中，K_Ⅰ＝0.56178，K_Ⅱ＝0.066490，则计算可得相当应力强度因子K_e＝0.56570<K_c＝0.575314，其中K_c为混凝土断裂韧度(它是用来描述混凝土自身断裂性能的参数，混凝土实际上是介于脆性与塑性之间的材料，随着混凝土强度的增加，其断裂韧度也会相应地降低，即表现为高强度混凝土比普通混凝土更“脆”)，因此第716次循环步及之后循环步为无效循环步，舍去。第715次循环步及之前循环步为有效循环步，并且有效循环步中的最后一步第715步的误差(K_e-K_c)/K_c＝1.484％，符合工程精度要求，故停止计算。此次计算即为所得的混凝土重力坝裂纹扩展路径，混凝土重力坝裂纹扩展路径图如图6所示，混凝土重力坝变形图如图7，其结果与工程经验吻合，即裂纹逐渐向混凝土重力坝坝体底部延伸。扩展完成后的混凝土重力坝裂纹尖端的x方向和y方向的局部应力场云图如图8所示，由图8其结果与断裂力学结论相吻合，即裂纹尖端处的应力集中现象明显并且裂纹尖端处的应力急剧增加。

本次计算所应用的混凝土重力坝模型结构简单，载荷及约束简化较多，所需要的计算成本也就较低，但计算效率及计算精度较高，并且由表可看出收敛速度较快。当混凝土重力坝的模型越接近实际工程模型时，模型结构也会相应地变得更加复杂，可以提高插值多项式的阶数来获得较高的精度和较快的收敛速度，最后得到满足工程精度要求的混凝土重力坝裂纹扩展路径。

由上表可以看出，本发明分析计算混凝土重力坝裂纹扩展路径时，前处理相对较少，收敛速度较快，而且通过较低的计算成本，获得了计算精度相对较高的混凝土重力坝裂纹扩展路径的数值解。

以上结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。

Claims (2)

1.一种基于高阶有限元法获取混凝土重力坝裂纹扩展路径的模拟计算方法，其特征在于：利用高阶有限元法计算出具有较高精度的混凝土重力坝应力场和位移场，进而得出高精度的应力强度因子及相当应力强度因子，最后获得相对准确的混凝土重力坝裂纹扩展路径，包括以下步骤：

步骤1、确定混凝土重力坝工程模型的几何尺寸，并创建混凝土重力坝的有限元模型；

步骤2、确定初始裂纹的位置及初始裂纹的长度，并将初始裂纹创建在混凝土重力坝有限元模型上；

步骤3、设置初始开裂方向、初始开裂位置及开裂步长；

步骤4、设置插值多项式阶数及循环次数；

每一次循环自动划分网格来进行混凝土重力坝裂纹扩展路径分析，得到每一循环步裂尖的应力强度因子，并求出相当应力强度因子，将相当应力强度因子和混凝土断裂韧度相比较，相当应力强度因子大于断裂韧度的循环步为有效循环，其余循环步舍去，判断有效循环步最后一步相当应力强度因子与混凝土断裂韧度之差是否在误差允许范围内，若不在误差允许范围内，则返回步骤3调小开裂步长并增大循环次数，继续进行计算；

所述步骤4具体包括以下几步：

步骤4.1、对插值多项式阶数和循环次数进行假定，其中插值多项式阶数为1～8；

步骤4.2、通过确定了插值多项式阶数的高阶有限元法计算得到含裂纹混凝土重力坝第一次循环的位移场和应力应变场，并且得到裂纹扩展角；

步骤4.3、利用围线积分法导出每一次循环的应力强度因子K_I和K_II；

步骤4.4、计算得到第一次循环的相当应力强度因子

应力强度因子K_I与K_II可以通过围线Γ₂，Γ₃及Γ₄上的线积分计算获得：

式中：

分别表示Ⅰ型与Ⅱ型应力强度因子的导出函数，

与

分别表示由

与

导出的牵引力向量；

步骤4.5、通过裂纹扩展角和裂纹扩展步长完成裂纹扩展第一步循环；

步骤4.6、开始第二步循环，重复步骤4.2至步骤4.6，完成裂纹扩展第二步循环，以此类推，直到完成所有循环次数，将所有计算得到的每一步循环的相当应力强度因子依次与混凝土断裂韧度对比；

若所有循环步的相当应力强度因子均大于混凝土断裂韧度对比，则增大循环次数，返回步骤4.1，若存在某些循环步的相当应力强度因子小于混凝土断裂韧度，则找出最接近且大于混凝土断裂韧度的那一步的相当应力强度因子，求得该步相当应力强度因子和混凝土断裂韧度的相对误差，看相对误差是否满足工程精度要求，如果满足，则裂纹扩展完成，如果不满足且相对误差与工程精度要求相差较小，则可尝试增大插值多项式阶数重复上述计算，如果不满足且相对误差与工程精度要求相差较大，则减小裂纹扩展步长并增大循环次数，减小后的裂纹扩展步长与增大后的循环次数的乘积应大于相当应力强度因子最接近且大于混凝土断裂韧度的那一步的循环步数与减小前的裂纹扩展步长的乘积，再重复上述步骤，直至满足工程精度要求。

2.根据权利要求1所述的基于高阶有限元法获取混凝土重力坝裂纹扩展路径的模拟计算方法，其特征在于：所述步骤1包括混凝土重力坝几何模型的建立，混凝土材料参数的设置，载荷及边界条件的施加。

Images