CN107832492B - 一种基于内聚力模型的钢结构腐蚀疲劳损伤计算方法 - Google Patents

Abstract

本发明属于钢结构腐蚀疲劳损伤分析技术领域，提供了一种基于内聚力模型的钢结构腐蚀疲劳损伤计算方法，步骤如下：定义一种内聚力单元模型；基于点蚀理论和蒙特卡洛方法处理腐蚀缺陷形成和扩展问题；利用本申请所定义的内聚力模型对拟研究结构进行时程分析；对拟研究结构模型网格进行动态更新；利用实时雨流计数法对实时应力进行处理，并计算等效应力幅和累积损伤量。通过本发明所述计算方法，可以分析腐蚀疲劳裂纹扩展过程，从而显著提高钢结构腐蚀疲劳损伤计算精度。

Description

一种基于内聚力模型的钢结构腐蚀疲劳损伤计算方法

技术领域

本发明属于钢结构腐蚀疲劳损伤分析技术领域，涉及一种利用有限元分析 钢结构腐蚀疲劳损伤的计算方法。

背景技术

腐蚀疲劳是钢结构在服役过程中不可避免出现的问题。工程中为简化计算， 通常采用Minner线性累积损伤准则来计算钢结构累积损伤：

其中，D表示累积损伤量；D_i表示第i种常幅应力循环造成的损伤，n_i表示 第i种常幅应力循环出现的次数，N_i表示该应力幅(S_i)对应的疲劳寿命。

应力幅值(S_i)可以通过经典雨流计数法得到。雨流计数法首先由M. Matsuishi和T.Endo提出。将应力-时间历程顺时针旋转90°，纵坐标表示时间， 横坐标表示应力，此时应力-时间历程和雨点从宝塔顶向下流动的情况相似。其 计数规则如下：

(1)以最高峰值或最低谷值(取绝对值较大者)为起点，重新安排应力-时间 历程。

(2)雨流依次从每个峰(谷)的内侧向下流，在下一个谷(峰)处落下，直 到对面有一个更低的谷值(或更高的峰值)停止。

(3)当雨流遇到来自上面屋顶流下的雨流时，即行停止。

(4)取出所有的全循环，并记录下各自的应力幅值(S_i)和均值。

应力幅(S_i)所对应的疲劳寿命(N_i)可以从试验数据回归得出的应力幅- 应力循环次数关系中得到。

由于受试验设备和试验环境的影响，试验钢结构构件尺寸受限，因此，该 方法无法预测任意尺寸钢结构构件的腐蚀疲劳寿命。

发明内容

本发明主要解决的技术问题是提供一种计算任意几何尺寸钢结构构件腐蚀 疲劳损伤的方法。该方法在考虑腐蚀、疲劳发生机理的基础上，分析钢结构构 件的腐蚀疲劳裂纹扩展过程，计算结构累积损伤，并估算其剩余疲劳寿命。

本发明的技术方案：

一种基于内聚力模型的钢结构腐蚀疲劳损伤计算方法，步骤如下：

(1)内聚力模型：内聚力单元模型如图2所示，内聚力单元节点对张开量以 上下界面

的位移变化||u||表示：||u||＝||u⁺-u^-||，其中u^±表示变形后上下界面的 位移；X_i(i＝1,2,...,n)表示三维空间的笛卡尔坐标，

表示变形后

对应点 的位置，

表示变形后中面Γ_d上对应点的位置；则有：

在中面Γ_d上一点P建立局部坐标系，以v_n,v_s,v_t表示局部坐标系在总体坐标 系中的方向余弦，构成正交旋转张量Θ_mi＝[v_n v_s v_t]，则单元节点在总体坐标系下 的位移变化向量：

Δ_m＝Θ_mi||u|| (3)

内聚力模型的本构关系以上下界面间粘结力τ_i和位移Δ_i在局部坐标系下之 间的函数关系表达：

τ_i＝τ(Δ_i) (4)

其张量表达式为：

其中，d∈[0,1]表示损伤变量，K为罚刚度，<-Δ₃>＝(Δ₃+|Δ₃|)/2；在混合加载模式下，内聚力模型的损伤准则以位移变化的形式表示：

r^t＝max{Δ⁰,maxλ^s},0≤s≤t (9)

d^t＝G(r^t) (10)

其中，λ表示界面间位移，λ^t、λ^s表示s时刻、t时刻对应的界面间位移，r^t表示t时刻损伤阈值。d^t表示单元损伤随时间的变化值，当d＝1时，单元破环， τ_i＝0。

在混合加载模式下使用能量释放率G来判断裂缝是否扩展，当能量释放率 G超过临界值G_c时界面开始分层，采用BK准则计算临界能量释放率：

G_c＝G_Ιc+(G_ΙΙc-G_Ιc)(G_shear/G)^η (11)

其中，η为材料参数，在混合加载模式下G＝G_Ι+G_shear,G_shear＝G_ΙΙ+G_ΙΙΙ。 G_Ι、G_ΙΙ、G_ΙΙΙ分别代表Ⅰ型、Ⅱ型、Ⅲ型断裂的能量释放率。G_Ιc、G_ΙΙc表示Ⅰ型、Ⅱ 型断裂的临界能量释放率。

(2)实时雨流计数法：实时雨流计数法主要用于提取变幅疲劳应力时程曲线 中的应力循环。通常情况下，应力时程曲线极为复杂，而雨流计数法分析需要 提供的应力必须为峰值与谷值相间的序列，因此，在处理真实情况下的应力序 列时必须先进行预处理，即删除应力序列中非峰谷值点。实时雨流计数法使用 两个动态堆栈，即应力峰值堆栈(datap)和应力谷值堆栈(datav)存放应力信息，按 照如下流程处理：

1)判断经预处理的应力堆栈中应力点数量是否大于1。若满足，则从该堆栈 中提取第一个应力点，置于峰值(谷值)堆栈内；否则，结束该峰值(谷值) 流程。

2)判断峰值(谷值)堆栈内应力点数量是否大于1。若满足，则将 Mxnew(Mnnew)赋值为峰值(谷值)堆栈内最后一个值，Mxold(Mnold)赋值为峰 值(谷值)堆栈内倒数第二个值；否则，进入步骤(3)。

3)重新判断预处理堆栈内应力点数量是否大于1，若满足，则跳转到谷值(峰 值)计算流程，否则，结束该峰值(谷值)计算流程。

4)判断是否满足Mxnew>Mxold(Mnnew>Mnold)。若满足，则进入步骤(5)； 否则，Mxold(Mnold)赋值为Mxnew(Mnnew)，进入步骤(3)。

5)判断谷值堆栈内应力点数量是否为1。若满足，则按照半循环处理，计算 应力幅、平均应力等应力循环信息，并统计应力点编号，从峰值(谷值)堆栈 移除应力点Mxold(Mnold)；否则按照全循环处理，计算应力幅、平均应力等 应力循环信息，并统计应力点编号，从峰值(谷值)堆栈移除应力点Mxold (Mnold)，从谷值堆栈移除应力点Mnnew(Mxnew)。跳转到步骤(3)。

(3)有限元模型内聚力单元更新：当裂缝前缘内聚力单元发生损伤时，查找 满足如下条件的实体单元：

1)与裂缝前缘内聚力单元外侧节点相连接；

2)实体单元的任意面不能与该内聚力单元任意界面重合；

3)满足条件步骤1)和2)的实体单元应当与另一个也满足条件步骤1)和 2)的实体单元存在公共面。

将满足以上条件且存在公共面的实体单元之间以内聚力单元连接。具体连 接方法为：假定存在公共面的两个实体单元分别为A和B。查询实体单元A和 B的所有节点编号及公共面上节点编号并保存；定义实体单元B在公共面上节 点集为C，B的其它节点组成节点集D，在节点集C的所有节点成员所在位置 新建重合节点，形成节点集E；删除实体单元B；以新建节点集E和节点集D 重新组合为新实体单元B；以节点集C和节点集E生成两个面，组成新的内聚 力单元。

(4)腐蚀效应处理：采用蒙特卡洛方法生成半球状点蚀蚀坑，随时间变化蚀 坑逐渐演化为半椭球状，应用有效应力集中系数K_f判定蚀坑发展的过程中新裂 纹的生成。

考虑蚀坑尺寸远小于构件尺寸，对于半椭球形蚀坑可以取应力集中系数 K_t＝1+2a/h。蚀坑底部曲率半径ρ＝a²/h，a为半椭球状蚀坑表面圆的半径，点 蚀蚀坑深度h随时间变化可以表示为：

h＝αt^β (13)

式中α、β为材料常数，二者可通过试验数据回归分析得到。

当K_f值为1～1.5时，点蚀蚀坑底部出现裂纹。新裂纹的处理参照单裂纹处理 方法。

本发明的有益效果：通过本发明所述计算方法，不但可以计算任意尺寸钢 结构腐蚀疲劳寿命，而且可以分析腐蚀疲劳裂纹扩展过程。

附图说明

图1是本发明方法实施流程图。

图2是内聚力单元模型示意图。

图3是实施例中实时雨流计数法的峰值堆栈处理方法框图。

图4是实施例中实时雨流计数法的谷值堆栈处理方法框图。

具体实施方式

下面结合附图对本发明的实施进行详细阐述，以使本发明的优点和特征能 更易于被本领域技术人员理解，从而对本发明的保护范围做出更为清楚明确的 界定。

本实施例应用Fortran语言对有限元软件进行二次开发，建立初始位置重合 的两个六节点三角形组成的12节点内聚力单元模型，结合实时雨流计数法、线 型累计损伤准则和采用蒙特卡洛方法生成的点蚀蚀坑及新裂纹，分析计算钢结 构腐蚀疲劳损伤。

内聚力单元的上下界面

上12个节点位移向量u_N为：

节点对之间的相对位移Δu_N为：

Δu_N＝Φ₁u_N (15)

其中，Φ₁＝(I_18×18|-I_18×18)，I_n×n为n×n单位阵，下同。

为了考虑节点对的形函数，引入自然坐标系，得到形函数表达式如下：

上下界面

上任意点(ξ，η)对应的界面相对位移为：

Δu(ξ,η)＝Φ₂Δu_N (17)

其中，Φ₂＝(φ₁I_3×3|φ₂I_3×3|φ₃I_3×3|φ₄I_3×3|φ₅I_3×3|φ₆I_3×3)。

假定节点初始坐标为

则上下界面节点对坐标为：

其中，Φ₃＝(I_18×18|I_18×18)，

从而，

根据中面坐标信息建立坐标转换阵如下：

本构关系和损伤准则参照“技术方案”章节进行处理。把内聚力单元模型 嵌入有限元分析软件计算分析钢结构应力。

使用实时雨流计数法分析从有限元软件接口提取的应力。根据“技术方案” 章节所述，在进行实时雨流计数法之前必须进行预处理。实时雨流计数法预处 理过程如下：取应力序列{σ_i}中的数据σ_i、σ_i+1和σ_i+2，当满足

(σ_i+1-σ_i)(σ_i+2-σ_i+1)＜0 (21)

时，点σ_i+1为峰值(谷值)点，σ_i、σ_i+1和σ_i+2均赋值为应力序列中下一个数据； 否则，删除数据点σ_i+1，σ_i+1和σ_i+2赋值为应力序列中下一个数据。重复上述操作， 直至σ_i+2赋值为应力序列中最后一个数据。

实时雨流计数法的实施流程参照“技术方案”章节、实时雨流计数法的峰值 堆栈处理方法框图(图3)和实时雨流计数法的谷值堆栈处理方法框图(图4) 进行。

每当可以采用实时雨流计数法提取出应力循环所对应的应力幅S_i时，通过此 时的内聚力单元面积之和∑A_CZE及按照裂缝前缘切平面方向计算出的剩余有效 面积预测值A_eff，计算累积损伤量：

则第i个应力循环造成的损伤量ΔD_i＝D_i-D_i-1，根据Minner线性累积损伤理论， 应力幅S_i所对应的疲劳寿命预测值为

随着S_i、N_i采集量的增大，可以得到S_i与N_i之间的回归关系，并用其预测剩 余腐蚀疲劳寿命。

完成实时雨流计数法分析后需要对有限元模型进行更新。有限元模型更新 包括两部分：(1)有限元模型内聚力单元更新；(2)有限元模型腐蚀蚀坑、裂 纹信息更新。

在每次进行有限元分析后，需要判断模型外表面处的内聚力单元是否形成 闭环。若是，则认为结构已经破环，计算结束；否则，进行实时雨流计数法分 析。完整的计算分析流程参见本发明方法实施流程图(图1)。

以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利 用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运 用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (1)

1.一种基于内聚力模型的钢结构腐蚀疲劳损伤计算方法，其特征在于，步骤如下：

(1)内聚力模型：内聚力单元节点对张开量以上下界面

的位移变化||u||表示：||u||＝||u⁺-u^-||，其中u^±表示变形后上、下界面的位移；X_i表示三维空间的笛卡尔坐标，i＝1,2,...,n；

表示变形后

对应点的位置，

表示变形后中面Γ_d上对应点的位置；则有：

在中面Γ_d上一点P建立局部坐标系，以v_n,v_s,v_t表示局部坐标系在总体坐标系中的方向余弦，构成正交旋转张量Θ_mi＝[v_n v_s v_t]，则单元节点在总体坐标系下的位移变化向量：

Δ_m＝Θ_mi||u|| (2)

内聚力模型的本构关系以上下界面间粘结力τ_i和位移Δ_i在局部坐标系下之间的函数关系表达：

τ_i＝τ(Δ_i) (3)

其张量表达式为：

其中，d∈[0,1]表示损伤变量，K为罚刚度，<-Δ₃>＝(Δ₃+|Δ₃|)/2；在混合加载模式下，内聚力模型的损伤准则以位移变化的形式表示：

r^t＝max{Δ⁰,maxλ^s},0≤s≤t (8)

d^t＝G(r^t) (9)

其中，λ表示界面间位移，λ^t、λ^s表示s时刻、t时刻对应的界面间位移，r^t表示t时刻损伤阈值；d^t表示单元损伤随时间的变化值，当d＝1时，单元破环，τ_i＝0；

在混合加载模式下，使用能量释放率G来判断裂缝是否扩展，当能量释放率G超过临界值G_c时界面开始分层，采用BK准则计算临界能量释放率：

G_c＝G_Ιc+(G_ΙΙc-G_Ιc)(G_shear/G)^η (10)

其中，η为材料参数，在混合加载模式下，G＝G_Ι+G_shear,G_shear＝G_ΙΙ+G_ΙΙΙ；G_Ι、G_ΙΙ、G_ΙΙΙ分别代表Ⅰ型、Ⅱ型、Ⅲ型断裂的能量释放率；G_Ιc、G_ΙΙc表示Ⅰ型、Ⅱ型断裂的临界能量释放率；

(2)实时雨流计数法：实时雨流计数法使用两个动态堆栈，即应力峰值堆栈和应力谷值堆栈存放应力信息，其中，应力峰值堆栈按照如下流程处理：

1)判断经预处理的应力堆栈中应力点数量是否大于1，若满足，则从应力堆栈中提取第一个应力点，置于应力峰值堆栈内；否则，结束应力峰值计算流程；

2)判断应力峰值堆栈内应力点数量是否大于1，若满足，则将Mxnew赋值为应力峰值堆栈内最后一个值，Mxold赋值为应力峰值堆栈内倒数第二个值；否则，进入步骤3)；

3)重新判断预处理的应力堆栈内应力点数量是否大于1，若满足，则跳转到应力谷值计算流程，否则，结束应力峰值计算流程；

4)判断是否满足Mxnew>Mxold，若满足，则进入步骤5)；否则，Mxold赋值为Mxnew，进入步骤3)；

5)判断应力谷值堆栈内应力点数量是否为1，若满足，则按照半循环处理，计算应力幅、平均应力的循环信息，并统计应力点编号，从应力峰值堆栈移除应力点Mxold；否则按照全循环处理，计算应力幅、平均应力的循环信息，并统计应力点编号，从应力峰值堆栈移除应力点Mxold，从应力谷值堆栈移除应力点Mnnew，跳转到步骤3)；

应力谷值堆栈按照如下流程处理：

1)判断经预处理的应力堆栈中应力点数量是否大于1，若满足，则从应力堆栈中提取第一个应力点，置于应力谷值堆栈内；否则，结束应力谷值计算流程；

2)判断应力谷值堆栈内应力点数量是否大于1，若满足，则将Mnnew赋值为应力谷值堆栈内最后一个值，Mnold赋值为应力谷值堆栈内倒数第二个值；否则，进入步骤3)；

3)重新判断预处理的应力堆栈内应力点数量是否大于1，若满足，则跳转到应力峰值计算流程，否则，结束应力谷值计算流程；

4)判断是否满足Mnnew>Mnold，若满足，则进入步骤5)；否则，Mnold赋值为Mnnew，进入步骤3)；

5)判断应力峰值堆栈内应力点数量是否为1，若满足，则按照半循环处理，计算应力幅、平均应力的循环信息，并统计应力点编号，从应力谷值堆栈移除应力点Mnold；否则按照全循环处理，计算应力幅、平均应力的循环信息，并统计应力点编号，从应力谷值堆栈移除应力点Mnold，从应力峰值堆栈移除应力点Mxnew，跳转到步骤3)；

(3)有限元模型内聚力单元更新：当裂缝前缘内聚力单元发生损伤时，查找满足如下条件的实体单元：

1)与裂缝前缘内聚力单元外侧节点相连接；

2)实体单元的任意面不能与该内聚力单元任意界面重合；

3)满足条件步骤1)和2)的实体单元应当与另一个也满足条件步骤1)和2)的实体单元存在公共面；

将满足以上条件且存在公共面的实体单元之间以内聚力单元连接，具体连接方法为：假定存在公共面的两个实体单元分别为A和B，查询实体单元A和B的所有节点编号及公共面上节点编号并保存；定义实体单元B在公共面上节点集为C，B的其它节点组成节点集D，在节点集C的所有节点成员所在位置新建重合节点，形成节点集E；删除实体单元B；以新建节点集E和节点集D重新组合为新实体单元B；以节点集C和节点集E生成两个面，组成新的内聚力单元；

(4)腐蚀效应处理：采用蒙特卡洛方法生成半球状点蚀蚀坑，随时间变化蚀坑逐渐演化为半椭球状，应用有效应力集中系数K_f判定蚀坑发展的过程中新裂纹的生成：

考虑蚀坑尺寸远小于构件尺寸，对于半椭球形蚀坑取应力集中系数K_t＝1+2a/h；蚀坑底部曲率半径ρ＝a²/h，a为半椭球状蚀坑表面圆的半径，点蚀蚀坑深度h随时间变化表示为：

h＝αt^β (12)

式中，α、β为材料常数，二者通过试验数据回归分析得到。

Images