CN111368482A - 循环荷载下钢筋与混凝土相互作用模拟计算方法 - Google Patents

Abstract

本发明公开了一种循环荷载下钢筋与混凝土相互作用模拟计算方法，建立包含钢筋与混凝土单元特性的三维有限元模型；给定预设的边界条件，施加循环往复荷载，确定计算时间步数与迭代步数；采用单弹簧联结单元法预测当前迭代步每个单元的位移值并计算钢筋与混凝土之间的相对滑移量；通过动态粘结滑移本构关系计算钢筋与混凝土之间的粘结强度；结合高斯点控制面积求解该节点处钢筋对混凝土的刚度贡献值，修正整体刚度矩阵；采用单弹簧联结单元法计算当前迭代步每个单元的位移与应力，更新钢筋与混凝土之间相对滑移量以及粘结强度。本发明准确计算循环荷载下钢筋与混凝土材料之间粘结强度，解决了现有方法很难模拟出钢筋的局部力和变形关系的问题。

Description

循环荷载下钢筋与混凝土相互作用模拟计算方法

技术领域

本发明涉及混凝土结构抗震分析计算，具体涉及一种循环荷载下钢筋与混凝土相互作用模拟计算方法。

背景技术

地震荷载作用下，高混凝土坝将不可避免的在坝体断面突变处及材料分区界面等部位出现易损区域，配筋已成为必不可少的工作，因此正确反映钢筋与混凝土之间的相互作用是仿真分析的重要方向。常规的钢筋混凝土有限元分析通常采用分布式模型、嵌入式模型或分离式模型，其中分布式模型假定钢筋均匀分布于混凝土实体单元内，通过实体单元积分点叠加钢筋刚度贡献的形式模拟钢筋混凝土相互作用；嵌入式模型则采用杆单元或梁单元模拟钢筋，通过位移约束条件强制钢筋单元与混凝土单元位移协同，最终将钢筋的刚度贡献等效叠加到混凝土节点处；分离式模型将分别离散钢筋单元与混凝土单元，在两者单元节点之间设置弹簧单元或接触单元，结合联结单元的粘结-滑移本构关系模拟钢筋混凝土之间的相互作用。

分布式模型和嵌入式模型假设钢筋混凝土不产生相对滑移，主观消除了钢筋位移自由度的运算，同时有限元模型中网格离散相对简单，计算效率能够得到提升，目前学者多采用此类模型进行水工结构抗震的仿真计算。但对于地震荷载而言，循环往复的荷载特性使得钢筋对混凝土的刚度贡献并不是一成不变的，显然分布式模型和嵌入式模型无法精准描述两者之间粘结滑移性能的变化过程，从而无法反映粘结性能变化对整体结构或局部构件承载能力和裂缝发展的影响。

发明内容

发明目的：本发明的目的是提供一种循环荷载下钢筋与混凝土相互作用模拟计算方法，清晰描述循环荷载下钢筋与混凝土材料之间粘结性能变化过程，解决了现有方法很难模拟出钢筋的局部力和变形关系的问题，并避免了人为选择法向刚度的困难。

技术方案：本发明所述的循环荷载下钢筋与混凝土相互作用模拟计算方法，包括：

(1)建立包含钢筋与混凝土单元特性的三维有限元模型，在钢筋单元节点处设置混凝土实体单元插值节点，并建立弹簧单元作为两者的联结单元；

(2)给定预设的边界条件，施加循环往复荷载，确定计算时间步数与迭代步数；

(3)进入当前运行步计算，采用单弹簧联结单元法预测当前迭代步每个单元的位移值并计算钢筋与混凝土之间的相对滑移量；

(4)根据步骤(3)得到的钢筋与混凝土之间的相对滑移量，通过动态粘结滑移本构关系计算钢筋与混凝土之间的粘结强度；

(5)根据步骤(4)得到的钢筋与混凝土之间的粘结强度，结合高斯点控制面积求解该节点处钢筋对混凝土的刚度贡献值，修正整体模型刚度矩阵；

(6)根据步骤(5)修正的整体刚度矩阵，采用单弹簧联结单元法计算当前迭代步每个单元的位移与应力，更新钢筋与混凝土之间相对滑移量以及粘结强度；

(7)根据步骤(6)得到的钢筋与混凝土之间相对滑移量以及粘结强度，结合高斯点控制面积更新钢筋对混凝土的刚度贡献值，更新节点荷载，计算出结构整体不平衡力；

(8)根据步骤(7)的计算结果，判断是否满足收敛条件，当未达到收敛条件时，重复步骤(5)—(8)计算下一迭代步；当达到收敛条件时，重复步骤(3)—(8)计算下一运行步，直至完所有运行步计算。

其中，所述步骤(3)中单弹簧联结单元法采用混合坐标系求解有限元方程，在局部坐标系中求解梁单元或杆单元的刚度矩阵，在整体坐标系中求解实体单元刚度矩阵，具体为：

根据公式

建立局部坐标系与整体坐标系之间的关系，其中：R为坐标系Oxyz与O'x'y'z'的坐标转换矩阵，

θ为绕Z轴的旋转角，Oxyz为混合坐标系中的整体坐标系，O'x'y'z'为梁单元局部坐标系，r为坐标系Oxyz与O^*x^*y^*z^*的坐标转换矩阵，O^*x^*y^*z^*梁单元节点局部坐标系用以协调相邻钢筋排列方向，通过节点相邻梁单元联合推求平均值，u，F分别为Oxyz中的位移向量和荷载向量，u'，F'分别为O'x'y'z'中的位移向量和荷载向量，u^*，F^*分别为O^*x^*y^*z^*中位移向量和荷载向量；

根据公式

计算局部坐标系下梁单元与节点单元的荷载向量，其中：k^*＝rkr^T＝rRk'Rr^T，k,k',k^*分别为Oxyz,O'x'y'z',O^*x^*y^*z^*中刚度矩阵；r^T表示坐标转换矩阵的转置矩阵；

根据公式

计算不协调网格中混凝土实体单元与钢筋单元交接点处钢筋节点的法向位移，其中：n为模型维数，m为钢筋插值节点周边的混凝土节点数，r_ij为插值系数，a为节点插值系数，u_j ^(k)为第k个插值点第j维的位移值；

根据增量形式的有限元方程

模拟钢筋与混凝土之间的相互作用，其中，f^*为局部坐标系O^*x^*y^*z^*中混凝土作用于钢筋的力，

f为粘结强度；s为粘结强度对应滑移量；D为钢筋直径；l为钢筋长度；d_s为当前滑移量；根据作用力的相互关系，-r^TΔf^*为整体坐标系Oxyz中钢筋作用于混凝土的力，Δu和ΔF表示Oxyz中的位移向量增量和荷载向量增量，Δu^*和ΔF^*表示O^*x^*y^*z^*中位移向量增量和荷载向量增量，Δf^*表示局部坐标系O^*x^*y^*z^*中混凝土作用于钢筋的力的增量。

所述步骤(4)中通过动态粘结-滑移本构关系计算粘结强度的过程具体如下：

节点处于单调加载阶段时，s_n≥s_n-1＝s_max，记录s_max＝s_n，根据模型公式

计算钢筋与混凝土之间滑移量与粘结强度的关系，其中：s_n为节点当前运行步滑移量，s_n-1为节点上一运行步滑移量，s_max为节点历史最大滑移量，τ为节点粘结强度，τ_s,τ_cr,τ_u,τ_r分别为胶着段，上升段，峰值段与退化段的粘结强度特征值，s为节点滑移量，s_s,s_cr,s_u,s_r分别为胶着段，上升段，峰值段与退化段的滑移量特征值，

k₁,k₂,k₃,k₄分别为胶着段，上升段，峰值段与退化段的斜率；

节点处于卸载阶段时，卸载路径分为两个部分：第一部分s_up≤s_n≤s_n-1≤s_max，根据公式τ_up＝-a_fτ_max计算残余粘结强度特征值，其中τ_up残余粘结强度；τ_max为节点历史最大粘结强度；a_f为粘结摩擦系数，取0.12，根据公式

计算残余滑移量特征点，此时粘结强度迅速下降至残余滑移量特征点，其中s_up为残余滑移量；s_max为节点历史最大滑移量；K_un为卸载刚度，取295N/m³，根据公式

计算节点粘结强度，第二部分s_n≤s_n-1≤s_up，此时粘结强度保持残余粘结强度τ_up至滑移零点。

所述步骤(4)中节点处于反向加卸载阶段时，粘结强度的计算与单调加卸载阶段的计算相同。

所述步骤(4)中节点处于多次循环后，根据公式

计算粘结强度退化系数，N为循环次数；节点处于稳定阶段s_n≥s_n-1≥s_r，粘结强度τ_n＝τ_n-1＝τ_r，由摩擦力提供并保持恒定。

所述步骤(6)中的单弹簧联结单元法和动态粘结-滑移本构关系计算过程与步骤(3)和步骤(4)中的计算过程相同。

有益效果：本发明准确描述了循环荷载下钢筋与混凝土材料之间粘结性能变化过程，解决了以往整体式钢筋混凝土动力模型很难模拟出钢筋的局部力和变形关系的困难，并且避免了人为选择法向刚度的困难，同时能够适应复杂的钢筋布置方案，在钢筋混凝土结构抗震分析中具有指导意义。

附图说明

图1是单弹簧联结单元法求解示意图；

图2是混合坐标系示意图；

图3是动态粘结-滑移本构关系骨架线示意图；

图4是动态粘结-滑移本构关系加载示意图；

图5是动态粘结-滑移本构关系卸载示意图；

图6是动态粘结-滑移本构关系反向加载示意图；

图7是动态粘结-滑移本构关系强度退化示意图；

图8是实施例一钢筋混凝土梁有限元示意图；

图9是实施例一循环荷载下侧表面钢筋节点粘结滑移关系示意图；

图10是实施例一循环荷载下钢筋0.5m处节点粘结滑移关系示意图；

图11是实施例一特征点受拉粘结峰值时钢筋变形沿长度分布示意图；

图12是实施例一特征点受拉粘结峰值时滑移量沿长度分布示意图；

图13是实施例一特征点受压粘结峰值时钢筋变形沿长度分布示意图；

图14是实施例一特征点受压粘结峰值时滑移量沿长度分布示意图；

具体实施方式

下面结合附图对发明进行进一步说明。

本发明公开的循环荷载下钢筋与混凝土相互作用模拟计算方法，包括以下步骤：

(1)建立包含钢筋与混凝土单元特性的三维有限元模型，如图1所示，在钢筋单元节点处设置混凝土实体单元插值节点，并建立弹簧单元作为两者的联结单元；

(2)给定合适的边界条件，施加合理的循环往复荷载，确定计算时间步数与迭代步数，进行数值计算；

(3)进入当前运行步计算，采用单弹簧联结单元法预测当前迭代步每个单元的位移值并计算钢筋与混凝土之间的相对滑移量，具体步骤如下：

3.1，建立局部坐标系与整体坐标系之间的关系，其公式为：

其中：R为坐标系Oxyz与O'x'y'z'的坐标转换矩阵，

θ为绕Z轴的旋转角；Oxyz为混合坐标系中的整体坐标系；O'x'y'z'为梁单元局部坐标系；r为坐标系Oxyz与O^*x^*y^*z^*的坐标转换矩阵；O^*x^*y^*z^*梁单元节点局部坐标系用以协调相邻钢筋排列方向，通过节点相邻梁单元联合推求平均值；u,F,u',F',u^*,F^*分别为Oxyz,O'x'y'z',O^*x^*y^*z^*中位移向量和荷载向量。

3.2，通过有限元方程计算局部坐标系下梁单元与节点单元的荷载向量，其公式为：

其中：k^*＝rkr^T＝rRk'Rr^T；k,k',k^*分别为Oxyz,O'x'y'z',O^*x^*y^*z^*中刚度矩阵，r^T表示坐标转换矩阵的转置矩阵；

3.3，通过插值的方式计算钢筋单元任意布置时不协调网格中混凝土实体单元与钢筋单元交接点处的插值节点位移，其公式为：

其中：m为混凝土单元被插值节点周边节点数。

3.4，通过插值的方式计算不协调网格中混凝土实体单元与钢筋单元交接点处钢筋节点的法向位移，其公式为：

其中：n为模型维数；r_ij为插值系数，即坐标转换矩阵元素；a为节点插值系数，即形状函数。

3.5，通过增量的形式建立有限元平衡方程模拟钢筋与混凝土之间的相互作用，其公式为：

其中：f^*为局部坐标系O^*x^*y^*z^*中混凝土作用于钢筋的力，

f为粘结强度，s为粘结强度对应滑移量；D为钢筋直径；l为钢筋长度；d_s为当前滑移量；根据作用力的相互关系，-r^TΔf^*为整体坐标系Oxyz中钢筋作用于混凝土的力。

(4)根据步骤(3)的计算结果，通过动态粘结滑移本构关系计算钢筋与混凝土之间的粘结强度，具体步骤如下：

4.1，当节点处于单调加载阶段时，即s_n≥s_n-1＝s_max，记录s_max＝s_n，根据《混凝土结构设计规范》(GB50010-2010)中的徐有邻五段线模型计算钢筋与混凝土之间滑移量与粘结强度的关系，其公式为：

其中：s_n为节点当前运行步滑移量；s_n-1为节点上一运行步滑移量；s_max为节点历史最大滑移量；τ,τ_s,τ_cr,τ_u,τ_r分别为节点粘结强度及各阶段粘结强度特征值；s,s_s,s_cr,s_u,s_r分别为节点滑移量及各阶段滑移量特征值；k₁,k₂,k₃,k₄分别为胶着段，上升段，峰值段与退化段的斜率；

4.2，当节点处于卸载阶段时，且s_up≤s_n≤s_n-1≤s_max，计算节点残余粘结强度特征值，其公式为：

τ_up＝-a_fτ_max

其中τ_up残余粘结强度；τ_max为节点历史最大粘结强度；a_f为粘结摩擦系数，取0.12。

4.3，此时粘结强度迅速下降至残余滑移量特征点，根据步骤4.2计算结果，计算节点残余滑移量特征点，其公式为：

其中s_up为残余滑移量；K_un为卸载刚度，根据腾智明试验经验参数取295N/m³。

4.4，根据步骤4.2和步骤4.3计算结果，计算节点粘结强度，其公式为：

4.5，当节点处于卸载阶段时，且s_n≤s_n-1≤s_up，粘结强度保持残余粘结强度τ_up至滑移零点。

4.6，当节点处于反向加卸载阶段时，因为钢筋与混凝土之间的粘结-滑移性能从机理上与施载方向无关，但外荷载的施加则有方向区分，因此其加卸载路径按原点中心对称。

4.7，当节点处于多次循环后，根据公式

计算粘结强度退化系数，其中：N为循环次数。

4.8，当节点处于稳定阶段时，即s_n≥s_n-1≥s_r，粘结强度τ_n＝τ_n-1＝τ_r，由摩擦力提供并保持恒定，与滑移量及加载方向无关。

(5)根据步骤(4)的计算结果，结合高斯点控制面积求解该节点处钢筋对混凝土的刚度贡献值，修正整体刚度矩阵；

(6)根据步骤(5)的计算结果，重复步骤3.1至步骤3.5计算当前迭代步每个单元的位移与应力，更新钢筋与混凝土之间相对滑移量，重复步骤4.1-4.8计算节点粘结强度；

(7)根据步骤(6)的计算结果，结合高斯点控制面积更新钢筋对混凝土的刚度贡献值，更新节点荷载，求解结构整体不平衡力；

(8)根据步骤(7)的计算结果，判断是否满足收敛条件。当未达到收敛条件时，重复步骤(5)至步骤(8)计算下一迭代步；当达到收敛条件时，重复步骤(3)至步骤(8)计算下一运行步，直至完所有运行步计算。

为了验证本方法的正确性，选取Shima的10号拉伸构件作为研究对象，圆柱体构件的三维有限元模型如图8所示，直径为0.5m，长度为0.76m，模型中心处设置一根直径为19mm的钢筋，在模型圆截面平行于坐标轴Y和坐标轴Z的外表面设置切向约束，钢筋两侧施加反方向的位移分级加载，并选取试件侧表面位置钢筋节点作为特征点以验证本文模型的正确性。

图9给出了特征点在循环荷载下的粘结滑移关系变化过程，可见所编译计算模块的仿真结果能够准确反映所嵌入的本构模型。

图10给出了循环荷载下钢筋0.5m处节点粘结滑移关系，由图可见循环荷载下构件中点附近钢筋与混凝土之间的粘结性能仅丧失了化学胶着力，机械咬合力仍发挥着决定作用，粘结强度仍处于上升阶段，与实际规律相符。

图11-14分别给出了特征点处于粘结强度峰值时钢筋变形与滑移量沿长度方向的分布情况，当端部荷载为拉加载时，钢筋各部位相对中点的轴向变形受粘结效应影响逐步降低，与实际规律相符。

图11给出了特征点受拉粘结峰值时钢筋变形与滑移量沿长度分布，由图可见粘结峰值时钢筋变形分布与Shima试验中一级加载与二级加载时钢筋变形分布形成了良好的过渡，与实际试验数据相契合。

通过图9-14可以看到，本发明提出一种循环荷载下钢筋与混凝土相互作用模拟计算方法，可以清晰描述循环荷载下钢筋与混凝土材料之间粘结性能变化过程，同时验证了本发明所提出分析方法的正确性和高效性。

Claims (6)

1.一种循环荷载下钢筋与混凝土相互作用模拟计算方法，其特征在于，包括：

(1)建立包含钢筋与混凝土单元特性的三维有限元模型，在钢筋单元节点处设置混凝土实体单元插值节点，并建立弹簧单元作为两者的联结单元；

(2)给定预设的边界条件，施加循环往复荷载，确定计算时间步数与迭代步数；

(3)进入当前运行步计算，采用单弹簧联结单元法预测当前迭代步每个单元的位移值并计算钢筋与混凝土之间的相对滑移量；

(4)根据步骤(3)得到的钢筋与混凝土之间的相对滑移量，通过动态粘结滑移本构关系计算钢筋与混凝土之间的粘结强度；

(5)根据步骤(4)得到的钢筋与混凝土之间的粘结强度，结合高斯点控制面积求解该节点处钢筋对混凝土的刚度贡献值，修正整体模型刚度矩阵；

(6)根据步骤(5)修正的整体刚度矩阵，采用单弹簧联结单元法计算当前迭代步每个单元的位移与应力，更新钢筋与混凝土之间相对滑移量以及粘结强度；

(7)根据步骤(6)得到的钢筋与混凝土之间相对滑移量以及粘结强度，结合高斯点控制面积更新钢筋对混凝土的刚度贡献值，更新节点荷载，计算出结构整体不平衡力；

(8)根据步骤(7)的计算结果，判断是否满足收敛条件，当未达到收敛条件时，重复步骤(5)—(8)计算下一迭代步；当达到收敛条件时，重复步骤(3)—(8)计算下一运行步，直至完所有运行步计算。

2.根据权利要求1所述的循环荷载下钢筋与混凝土相互作用模拟计算方法，其特征在于，所述步骤(3)中单弹簧联结单元法采用混合坐标系求解有限元方程，在局部坐标系中求解梁单元或杆单元的刚度矩阵，在整体坐标系中求解实体单元刚度矩阵，具体为：

根据公式

建立局部坐标系与整体坐标系之间的关系，其中：R为坐标系Oxyz与O'x'y'z'的坐标转换矩阵，

θ为绕Z轴的旋转角，Oxyz为混合坐标系中的整体坐标系，O'x'y'z'为梁单元局部坐标系，r为坐标系Oxyz与O^*x^*y^*z^*的坐标转换矩阵，O^*x^*y^*z^*梁单元节点局部坐标系用以协调相邻钢筋排列方向，通过节点相邻梁单元联合推求平均值，u，F分别为Oxyz中的位移向量和荷载向量，u'，F'分别为O'x'y'z'中的位移向量和荷载向量，u^*，F^*分别为O^*x^*y^*z^*中位移向量和荷载向量；

根据公式

计算局部坐标系下梁单元与节点单元的荷载向量，其中：k^*＝rkr^T＝rRk'Rr^T，k,k',k^*分别为Oxyz,O'x'y'z',O^*x^*y^*z^*中刚度矩阵；r^T表示坐标转换矩阵的转置矩阵；

根据公式

计算不协调网格中混凝土实体单元与钢筋单元交接点处钢筋节点的法向位移，其中：n为模型维数，m为钢筋插值节点周边的混凝土节点数，r_ij为插值系数，a为节点插值系数，u_j ^(k)为第k个插值点第j维的位移值；

根据增量形式的有限元方程

模拟钢筋与混凝土之间的相互作用，其中，f^*为局部坐标系O^*x^*y^*z中混凝土作用于钢筋的力，

f为粘结强度；s为粘结强度对应滑移量；D为钢筋直径；l为钢筋长度；d_s为当前滑移量；根据作用力的相互关系，-r^TΔf^*为整体坐标系Oxyz中钢筋作用于混凝土的力，Δu和ΔF表示Oxyz中的位移向量增量和荷载向量增量，Δu^*和ΔF^*表示O^*x^*y^*z^*中位移向量增量和荷载向量增量，Δf^*表示局部坐标系O^*x^*y^*z^*中混凝土作用于钢筋的力的增量。

3.根据权利要求1所述的循环荷载下钢筋与混凝土相互作用模拟计算方法，其特征在于，所述步骤(4)中通过动态粘结-滑移本构关系计算粘结强度的过程具体如下：

节点处于单调加载阶段时，s_n≥s_n-1＝s_max，记录s_max＝s_n，根据模型公式

计算钢筋与混凝土之间滑移量与粘结强度的关系，其中：s_n为节点当前运行步滑移量，s_n-1为节点上一运行步滑移量，s_max为节点历史最大滑移量，τ为节点粘结强度，τ_s,τ_cr,τ_u,τ_r分别为胶着段，上升段，峰值段与退化段的粘结强度特征值，s为节点滑移量，s_s,s_cr,s_u,s_r分别为胶着段，上升段，峰值段与退化段的滑移量特征值，

k₁,k₂,k₃,k₄分别为胶着段，上升段，峰值段与退化段的斜率；

节点处于卸载阶段时，卸载路径分为两个部分：第一部分s_up≤s_n≤s_n-1≤s_max，根据公式τ_up＝-a_fτ_max计算残余粘结强度特征值，其中τ_up残余粘结强度；τ_max为节点历史最大粘结强度；a_f为粘结摩擦系数，取0.12，根据公式

计算残余滑移量特征点，此时粘结强度迅速下降至残余滑移量特征点，其中s_up为残余滑移量；s_max为节点历史最大滑移量；K_un为卸载刚度，取295N/m³，根据公式

计算节点粘结强度，第二部分s_n≤s_n-1≤s_up，此时粘结强度保持残余粘结强度τ_up至滑移零点。

4.根据权利要求3所述的循环荷载下钢筋与混凝土相互作用模拟计算方法，其特征在于，所述步骤(4)中节点处于反向加卸载阶段时，粘结强度的计算与单调加卸载阶段的计算相同。

5.根据权利要求3所述的循环荷载下钢筋与混凝土相互作用模拟计算方法，其特征在于，所述步骤(4)中节点处于多次循环后，根据公式

计算粘结强度退化系数，N为循环次数；节点处于稳定阶段s_n≥s_n-1≥s_r，粘结强度τ_n＝τ_n-1＝τ_r，由摩擦力提供并保持恒定。

6.根据权利要求1所述的循环荷载下钢筋与混凝土相互作用模拟计算方法，其特征在于，所述步骤(6)中的单弹簧联结单元法和动态粘结-滑移本构关系计算过程与步骤(3)和步骤(4)中的计算过程相同。

Images