裂纹尖端有限元模型建模技术研究和精度分析
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
裂纹尖端有限元模型建模技术研究和精度分析
王予津 1，张俊芳 2 1.上海交通大学机械与动力工程学院，200240 2.上海交通大学材料科学与工程学院，200240 [ 摘 要 ] 在进行有限元模拟裂纹扩展过程中，需要建立裂纹尖端模型、提取裂纹应力强度 因子等参数和裂纹扩展速率和方向的计算。 如果模拟某一部件从裂纹萌生到部件断裂这一过程， 就需要重复对不同时期的裂纹进行建模和分析。 本文基于于 ANSYS/APDL 语言， 通过建立模块 化的裂纹尖端模型，从而实现裂纹尖端模型的快速重建。并在上述技术基础上，又采用位移外 推法，J 积分和虚拟裂纹闭合法计算裂纹的应力强度因子，并与解析解进行了对比，分析不同裂 纹深度下，各个计算方法的误差。
1
前言
在进行断裂力学分析中，最大的困难是模拟裂纹从萌生到断裂的整个过程。这是在一个疲 劳载荷周期内，裂纹的扩展速率、开裂角度以及支裂纹等都会影响下一阶段的断裂模型的建立。 本文基于 ANSYS/APDL 语言提出了一种裂纹尖端模块化建模方法， 将已建立好的裂纹模块， 通 过各种形式的组和，快速得到所需裂纹形式。
K1_K1 K1_K2 K1_K3 K1_J_average K1_J_average k1_ext
位移外推法，图 9（a） 位移外推法，图 9（b） 位移外推法，图 9（c） J 积分（闭合面中间位置半径作为积 分路径） J 积分（闭合面外沿半径作为积分路 径） 虚拟裂纹闭合法
对平板不同裂纹深度进行分析，可以得到解析解的应力数值，同时也给出了各个数值计算 结果与解析解之间的相对误差，如图 10 所示，表 3 给出的数值分析结果与解析解之间的误差， 其中，误差分析依据的解析解为式 7 所给的结果。
图 10 不同裂纹深度下几种算法所得应力强度因子 表 3 数值分析结果与解析解之间的误差
2W
算例中，基于有限元分析方法，采用位移外推法、J 积分和虚拟裂纹闭合法。其中位移外推 法选取三个不同形式的节点进行计算。如图 9 所示，箭头所指为选取计算的节点。
a.三个节点选取一个单元内
b.三个节点选取相邻两个单元内 图 9 位移外推法中节点的选取
c.三个节点选取裂纹面上
J 积分选取图 9(c)中除尖端节点以外两个节点半径作为积分路径的单元积分路径。 虚拟裂纹闭合法选取的节点为裂纹尖端处裂纹闭合面上的两个奇异单元对应的角节点（非 尖端节点）。 3.4 求解精度分析 算例中，将各种方法和求解结果进行对照变化，其关系如表 2 所示。 表 2 结果与算法名称对照表 计算方法 K1_value1 K1_value2 K1_value3 对应名称 解析解，如式 5 解析解，如式 6 解析解，如式 7
（2）
图 2 位移外推法采用的裂纹尖端有限元模型 在使用位移外推法时，需要设定裂纹尖端处的节点和 A，B 以及，裂纹面 1.2 J 积分； J 积分是从能量守恒的角度去考虑裂纹在开裂过程中所释放的能量， 即所释放的能量应当为 围道上单元的应变能-外力作用在边界上的功。因此，在进行围道积分过程中，无需考虑应力的 奇异，但应当尽量选取一个式中的积分路径，即不能离裂纹尖端处过近，也不能太远。其 J 积 分的计算形式如式所示 3 所示，图 3 是式 3 对应的积分路径。
图 5 分析所需的各种裂纹模块 2.2 裂纹扩展模拟 裂纹扩展模拟的具体分析流程如图 2 所示。将含有裂纹的模型建立之后，加载边界条件， 进行静力学分析， 得出当前模型下的裂纹应力强度因子； 然后， 将应力强度因子的数值带入 Paris 公式中计算裂纹的扩展速率，并换算出下一阶段裂纹的扩展距离；同时通过对裂纹尖端处单元 的应变能进行评估，判断单元应变能最大值所在的方向，将其作为裂纹开裂的方向。最终判断 裂纹扩展距离 a 是否达到临界应力强度因子所对应的裂纹深度，从而确定模型是否发生断裂。 通过上述流程， 就可以很快的通过 ANSYS/APDL 语言进行所需的裂纹参数计算， 并且考虑 裂纹扩展速率和开裂角等参数，可以快速建立新的裂纹模型，从而实现裂纹扩展的预测。
（4）
图 4 虚拟裂纹闭合法
2
裂纹模型模块化
2.1 裂纹模块的建立
当前，断裂力学分析软件主要有基于有限元法的 ZENCRACK 和基于边界元法的 BEASY。 前者主要的思路是采用 crack block 技术，将部件的有限元模型上可能存在裂纹的部位，用其 crack block 进行替换， 从而实现模型的建立。 后者， 目前只能进行二维的断裂力学。 当前， Ansys 也能够进行断裂力学应力强度因子，J 积分，能量释放率的计算，但是该方法对于复杂结构的裂 纹建立却是很难实现。 本文采用了与 Zencrack 相同的思路，将可能存在的裂纹形式进行模块化建立，针对不同的 裂纹进行模块组和和替换，并对裂纹模块中的指定节点进行数据提取，用于后续的裂纹应力强 度因子和 J 积分的计算。 首先， 建立裂纹尖端的模块化数据库， 该数据库采用 ANSYS/APDL 语言建立， 如图 5 所示。 其次，建立无裂纹模型。然后，指定部件存在裂纹的部位，将裂纹模块替换正常的六面体单元， 保证裂纹模块的角点和单元的角点坐标重合，编号压缩。最终实现裂纹尖端模型的建立。
同时，断裂力学中裂纹应力强度因子等参数的计算常采用的方法有限元法、边界元法以及 基于哈密顿体系的辛几何算法。前者属于数值解法，在进行裂纹应力强度因子计算时，存在应 力奇异的问题；而后两种方法属于半解析半数值的方法，能够很好的解决有限元法在裂纹尖端 存在应力奇异问题。但是，由于裂纹尖端常常有一定范围的塑性区域，该部位属于材料非线性 问题，边界元法和辛几何算法是无法考虑求解域内的材料塑性问题，这也是边界元法和哈密顿 体系下的力学分析无法在工程中广泛应用的一个重要原因之一。 基于有限元法的断裂力学分析可分为直接法和间接法。直接法指通过裂纹尖端附近的应力 场和位移场的关系直接确定应力强度因子的值；间接法指通过有限元法分析，依据能量守恒的 原理，求出中间变量（如：J 积分），进而根据应力强度因子与这些中间变量之间的关系，简介 换算出应力强度因子的数值。 直接法有自身的一些劣势，式 1 所列的是直裂纹尖端的应力场分布公式。KI 为应力强度因 子，KI 为一类裂纹的应力强度因子，r 是观测点与裂纹尖端之间的距离，可见随着观察点与裂 纹尖端节点距离的减小，应力值会随之增大，并趋近于无穷大。
[ 关键词 ] 裂纹扩展预测，有限元法，裂纹模型重建，精度分析
Crack Tip Modeling Technology Research and Precision Analysis
Yujin Wang, Junfang Zhang 1.School of Mechanical Engineering,Shanghai Jiaotong University,200240 2.School of Material Science and Engineering,Shanghai Jiaotong University,200240 [ Abstract ] There are three steps in the simulation of crack propagation process with FEM, the first step is establishing the model of the crack tip , the second step is extracting the stress intensity factor and other parameters, the last step is calculating the crack growth rate and propagation direction. If we want to simulate a part from crack initiation to fracture, we need to repeat the simulation with the different periods of crack. Based on ANSYS / APDL language, through the establishment of a modular crack tip model, in order to achieve the rapid reconstruction of the crack tip model. And on the basis of the literal, but also by the displacement extrapolation, J integral and virtual crack closure calculation of the stress intensity factor, and whose solution were compared with the analytical solution. [ Keyword ] Crack propagation prediction, Finite element method, Crack model reconstruction, accuracy analysis
J dy T
u ds x
（3）
i j
其中，
0
ij
d ij 。
图 3 J 积分路径 1.3 虚拟裂纹闭合法
虚拟裂纹闭合法也是从能量的角度去分析裂纹的扩展，假设裂纹在下一个阶段扩展距离为 裂纹尖端处单元的长度，分析其开裂所需要的能量。 相对于围道积分而言，其形式较为简单，如式 4 所示，图 4 是其裂纹尖端的模型。其中 Fy1 是沿着 y 向的外部载荷，Δ v3，4 是 3 和 4 节点的变形后的距离，B 是板的厚度，Δ a 为虚拟的裂 纹扩展距离，一般设为裂纹尖端处单元的长度。
图 6 裂纹扩展模拟流程 分析中尖端附近的单元采用 1/4 节点法的奇异单元进行模型建立，并且围绕尖端部位的单 元层数和周向的单元个数可以进行人为控制，从而可以保证分析所得结果的精度。
3
算例
平板中心裂纹模型如图 7 所示，具体几何尺寸如表 1 所示，承受垂直方向的拉应力 0.5641895psi。其有限元模型选取为平板裂纹的 1/4 模型，如图 8 所示。
KI Fra Baidu bibliotek a
2W a tan a 2W
（5）
2 4 a a K I a sec 1 0.025 0.06 2W W W
a a
（6） （7）
K I a sec
3.3 有限元法
图 1 1/4 节点法下的奇异单元中间节点位置 1.1 位移外推法 位移外推法所采用的计算形式如式 2 所示。图 2 是该公式对应的裂纹尖端处有限元模型。 其中 uA,B 为 1/4 节点和单元角节点位移差值，r 为 1/4 节点到裂纹尖端节点的距离。
KI
2 E u A, B 4(1 v 2 ) r
KI 2 r
（1）
1975 年，Hanshell 提出了 1/4 节点法，就是将单元中间节点移动到 1/4 节点位置，该单元的 插值形函数也会发生变化，会使得中间节点的位移具有奇异性。因此，在进行计算时，中间节 点应力的奇异性和它的位移奇异性正好抵消，而不会出现应力奇异[1]，这一方法在断裂力学直 接法——位移外推法得到了应用。奇异单元如图 1 所示
裂纹
图 7 平板中心裂纹模型
图 8 1/4 平板裂纹模型 表 1 平板中心裂纹几何参数 几何参数 A B H T EX 数值 1in 5in 5in 0.25in 3e7psi
PRXY
0.3
采用三种位移外推法、 J 积分和虚拟裂纹闭合法和解析解对不同裂纹深度模型进行应力强度 因子计算误差的分析。 3.2 裂纹尖端应力强度因子解析解 本文选用了三应力强度因子解析解模型。如式 5~7 所示[2-4]。W 为平板宽度， a 为裂纹深度。
王予津 1，张俊芳 2 1.上海交通大学机械与动力工程学院，200240 2.上海交通大学材料科学与工程学院，200240 [ 摘 要 ] 在进行有限元模拟裂纹扩展过程中，需要建立裂纹尖端模型、提取裂纹应力强度 因子等参数和裂纹扩展速率和方向的计算。 如果模拟某一部件从裂纹萌生到部件断裂这一过程， 就需要重复对不同时期的裂纹进行建模和分析。 本文基于于 ANSYS/APDL 语言， 通过建立模块 化的裂纹尖端模型，从而实现裂纹尖端模型的快速重建。并在上述技术基础上，又采用位移外 推法，J 积分和虚拟裂纹闭合法计算裂纹的应力强度因子，并与解析解进行了对比，分析不同裂 纹深度下，各个计算方法的误差。
1
前言
在进行断裂力学分析中，最大的困难是模拟裂纹从萌生到断裂的整个过程。这是在一个疲 劳载荷周期内，裂纹的扩展速率、开裂角度以及支裂纹等都会影响下一阶段的断裂模型的建立。 本文基于 ANSYS/APDL 语言提出了一种裂纹尖端模块化建模方法， 将已建立好的裂纹模块， 通 过各种形式的组和，快速得到所需裂纹形式。
K1_K1 K1_K2 K1_K3 K1_J_average K1_J_average k1_ext
位移外推法，图 9（a） 位移外推法，图 9（b） 位移外推法，图 9（c） J 积分（闭合面中间位置半径作为积 分路径） J 积分（闭合面外沿半径作为积分路 径） 虚拟裂纹闭合法
对平板不同裂纹深度进行分析，可以得到解析解的应力数值，同时也给出了各个数值计算 结果与解析解之间的相对误差，如图 10 所示，表 3 给出的数值分析结果与解析解之间的误差， 其中，误差分析依据的解析解为式 7 所给的结果。
图 10 不同裂纹深度下几种算法所得应力强度因子 表 3 数值分析结果与解析解之间的误差
2W
算例中，基于有限元分析方法，采用位移外推法、J 积分和虚拟裂纹闭合法。其中位移外推 法选取三个不同形式的节点进行计算。如图 9 所示，箭头所指为选取计算的节点。
a.三个节点选取一个单元内
b.三个节点选取相邻两个单元内 图 9 位移外推法中节点的选取
c.三个节点选取裂纹面上
J 积分选取图 9(c)中除尖端节点以外两个节点半径作为积分路径的单元积分路径。 虚拟裂纹闭合法选取的节点为裂纹尖端处裂纹闭合面上的两个奇异单元对应的角节点（非 尖端节点）。 3.4 求解精度分析 算例中，将各种方法和求解结果进行对照变化，其关系如表 2 所示。 表 2 结果与算法名称对照表 计算方法 K1_value1 K1_value2 K1_value3 对应名称 解析解，如式 5 解析解，如式 6 解析解，如式 7
（2）
图 2 位移外推法采用的裂纹尖端有限元模型 在使用位移外推法时，需要设定裂纹尖端处的节点和 A，B 以及，裂纹面 1.2 J 积分； J 积分是从能量守恒的角度去考虑裂纹在开裂过程中所释放的能量， 即所释放的能量应当为 围道上单元的应变能-外力作用在边界上的功。因此，在进行围道积分过程中，无需考虑应力的 奇异，但应当尽量选取一个式中的积分路径，即不能离裂纹尖端处过近，也不能太远。其 J 积 分的计算形式如式所示 3 所示，图 3 是式 3 对应的积分路径。
图 5 分析所需的各种裂纹模块 2.2 裂纹扩展模拟 裂纹扩展模拟的具体分析流程如图 2 所示。将含有裂纹的模型建立之后，加载边界条件， 进行静力学分析， 得出当前模型下的裂纹应力强度因子； 然后， 将应力强度因子的数值带入 Paris 公式中计算裂纹的扩展速率，并换算出下一阶段裂纹的扩展距离；同时通过对裂纹尖端处单元 的应变能进行评估，判断单元应变能最大值所在的方向，将其作为裂纹开裂的方向。最终判断 裂纹扩展距离 a 是否达到临界应力强度因子所对应的裂纹深度，从而确定模型是否发生断裂。 通过上述流程， 就可以很快的通过 ANSYS/APDL 语言进行所需的裂纹参数计算， 并且考虑 裂纹扩展速率和开裂角等参数，可以快速建立新的裂纹模型，从而实现裂纹扩展的预测。
（4）
图 4 虚拟裂纹闭合法
2
裂纹模型模块化
2.1 裂纹模块的建立
当前，断裂力学分析软件主要有基于有限元法的 ZENCRACK 和基于边界元法的 BEASY。 前者主要的思路是采用 crack block 技术，将部件的有限元模型上可能存在裂纹的部位，用其 crack block 进行替换， 从而实现模型的建立。 后者， 目前只能进行二维的断裂力学。 当前， Ansys 也能够进行断裂力学应力强度因子，J 积分，能量释放率的计算，但是该方法对于复杂结构的裂 纹建立却是很难实现。 本文采用了与 Zencrack 相同的思路，将可能存在的裂纹形式进行模块化建立，针对不同的 裂纹进行模块组和和替换，并对裂纹模块中的指定节点进行数据提取，用于后续的裂纹应力强 度因子和 J 积分的计算。 首先， 建立裂纹尖端的模块化数据库， 该数据库采用 ANSYS/APDL 语言建立， 如图 5 所示。 其次，建立无裂纹模型。然后，指定部件存在裂纹的部位，将裂纹模块替换正常的六面体单元， 保证裂纹模块的角点和单元的角点坐标重合，编号压缩。最终实现裂纹尖端模型的建立。
同时，断裂力学中裂纹应力强度因子等参数的计算常采用的方法有限元法、边界元法以及 基于哈密顿体系的辛几何算法。前者属于数值解法，在进行裂纹应力强度因子计算时，存在应 力奇异的问题；而后两种方法属于半解析半数值的方法，能够很好的解决有限元法在裂纹尖端 存在应力奇异问题。但是，由于裂纹尖端常常有一定范围的塑性区域，该部位属于材料非线性 问题，边界元法和辛几何算法是无法考虑求解域内的材料塑性问题，这也是边界元法和哈密顿 体系下的力学分析无法在工程中广泛应用的一个重要原因之一。 基于有限元法的断裂力学分析可分为直接法和间接法。直接法指通过裂纹尖端附近的应力 场和位移场的关系直接确定应力强度因子的值；间接法指通过有限元法分析，依据能量守恒的 原理，求出中间变量（如：J 积分），进而根据应力强度因子与这些中间变量之间的关系，简介 换算出应力强度因子的数值。 直接法有自身的一些劣势，式 1 所列的是直裂纹尖端的应力场分布公式。KI 为应力强度因 子，KI 为一类裂纹的应力强度因子，r 是观测点与裂纹尖端之间的距离，可见随着观察点与裂 纹尖端节点距离的减小，应力值会随之增大，并趋近于无穷大。
[ 关键词 ] 裂纹扩展预测，有限元法，裂纹模型重建，精度分析
Crack Tip Modeling Technology Research and Precision Analysis
Yujin Wang, Junfang Zhang 1.School of Mechanical Engineering,Shanghai Jiaotong University,200240 2.School of Material Science and Engineering,Shanghai Jiaotong University,200240 [ Abstract ] There are three steps in the simulation of crack propagation process with FEM, the first step is establishing the model of the crack tip , the second step is extracting the stress intensity factor and other parameters, the last step is calculating the crack growth rate and propagation direction. If we want to simulate a part from crack initiation to fracture, we need to repeat the simulation with the different periods of crack. Based on ANSYS / APDL language, through the establishment of a modular crack tip model, in order to achieve the rapid reconstruction of the crack tip model. And on the basis of the literal, but also by the displacement extrapolation, J integral and virtual crack closure calculation of the stress intensity factor, and whose solution were compared with the analytical solution. [ Keyword ] Crack propagation prediction, Finite element method, Crack model reconstruction, accuracy analysis
J dy T
u ds x
（3）
i j
其中，
0
ij
d ij 。
图 3 J 积分路径 1.3 虚拟裂纹闭合法
虚拟裂纹闭合法也是从能量的角度去分析裂纹的扩展，假设裂纹在下一个阶段扩展距离为 裂纹尖端处单元的长度，分析其开裂所需要的能量。 相对于围道积分而言，其形式较为简单，如式 4 所示，图 4 是其裂纹尖端的模型。其中 Fy1 是沿着 y 向的外部载荷，Δ v3，4 是 3 和 4 节点的变形后的距离，B 是板的厚度，Δ a 为虚拟的裂 纹扩展距离，一般设为裂纹尖端处单元的长度。
图 6 裂纹扩展模拟流程 分析中尖端附近的单元采用 1/4 节点法的奇异单元进行模型建立，并且围绕尖端部位的单 元层数和周向的单元个数可以进行人为控制，从而可以保证分析所得结果的精度。
3
算例
平板中心裂纹模型如图 7 所示，具体几何尺寸如表 1 所示，承受垂直方向的拉应力 0.5641895psi。其有限元模型选取为平板裂纹的 1/4 模型，如图 8 所示。
KI Fra Baidu bibliotek a
2W a tan a 2W
（5）
2 4 a a K I a sec 1 0.025 0.06 2W W W
a a
（6） （7）
K I a sec
3.3 有限元法
图 1 1/4 节点法下的奇异单元中间节点位置 1.1 位移外推法 位移外推法所采用的计算形式如式 2 所示。图 2 是该公式对应的裂纹尖端处有限元模型。 其中 uA,B 为 1/4 节点和单元角节点位移差值，r 为 1/4 节点到裂纹尖端节点的距离。
KI
2 E u A, B 4(1 v 2 ) r
KI 2 r
（1）
1975 年，Hanshell 提出了 1/4 节点法，就是将单元中间节点移动到 1/4 节点位置，该单元的 插值形函数也会发生变化，会使得中间节点的位移具有奇异性。因此，在进行计算时，中间节 点应力的奇异性和它的位移奇异性正好抵消，而不会出现应力奇异[1]，这一方法在断裂力学直 接法——位移外推法得到了应用。奇异单元如图 1 所示
裂纹
图 7 平板中心裂纹模型
图 8 1/4 平板裂纹模型 表 1 平板中心裂纹几何参数 几何参数 A B H T EX 数值 1in 5in 5in 0.25in 3e7psi
PRXY
0.3
采用三种位移外推法、 J 积分和虚拟裂纹闭合法和解析解对不同裂纹深度模型进行应力强度 因子计算误差的分析。 3.2 裂纹尖端应力强度因子解析解 本文选用了三应力强度因子解析解模型。如式 5~7 所示[2-4]。W 为平板宽度， a 为裂纹深度。