基于Abaqus求解考虑裂纹面接触的裂尖应力强度因子的可行性分析
对于承受滚动接触载荷的含裂纹结构，如车轮、轴承等，由于压缩载荷的作用，裂纹面有可能会出现闭合的情况，为了准确考虑裂纹面闭合对于疲劳裂纹扩展的影响，往往需要在考虑裂纹面接触的情况下求解裂纹尖端的应力强度因子。
通用有限元软件Abaqus具备求解含裂纹结构在任意载荷作用下的应力强度因子K、J积分和T应力等断裂参量的功能。然而，一些研究表明，Abaqus无法准确求解考虑裂纹面接触时裂纹尖端的应力强度因子。为此，本文以承受滚动接触载荷作用的含斜裂纹平板为例，分别基于有限元软件Abaqus、断裂分析软件Franc3D在考虑和不考虑裂纹面接触的情况下求解裂纹尖端的应力强度因子，并分析了在此类情况下采用Abaqus准确求解裂尖应力强度因子的可行性。
有限元建模
有限元建模
模型简介
考虑长度为30 mm，高度和厚度均为10 mm的平板，在平板中部存在一条长度为4 mm的穿透型斜裂纹，如图1所示。假设平板为线弹性材料，材料的弹性模量E为200 GPa，泊松比ν为0.3。
图1 含斜裂纹平板示意图
该平板受到滚动接触载荷的作用，当滚动接触载荷移动到斜裂纹上方所在的区域时，裂纹面将会产生接触。在本例中，将该滚动接触载荷简化为Hertz接触。因此，作用在平板上的接触压力p(x)可表示为：
式中，p0为接触中心的最大接触压力，本例中取为100 MPa；a为接触斑半径，取为0.5mm。载荷的移动速度v取为30 mm/s。
下面，需要计算在该滚动接触载荷作用下，裂纹尖端应力强度因子的变化。
基于有限元软件Abaqus的建模方法
当采用有限元软件Abaqus进行建模时，需建立含有裂纹的平板有限元模型，并通过内置的围线积分法求解裂纹体的应力强度因子。
为保证计算精度，在平板表面和裂纹尖端区域采用细密网格进行划分，网格尺寸分别取为0.2 mm和0.1 mm，对于远离平板表面和裂纹尖端的区域，则采用较为稀疏的网格划分。
为模拟裂纹面间的接触行为，本文在上、下裂纹面之间建立自接触（Self contact），并取裂纹面间的摩擦系数为0.1。约束平板底面的全部位移自由度，并通过用户子程序DLOAD施加滚动接触载荷。最终建立的有限元模型如图2所示。
图2 含斜裂纹平板示意图（有限元软件Abaqus）
基于断裂分析软件Franc3D的建模方法
当采用断裂分析软件Franc3D求解应力强度因子时，需首先在有限元软件Abaqus中建立不含裂纹的平板模型，并设置好载荷和边界条件。然后，将不含裂纹的有限元模型导入到Franc3D中，并在Franc3D中植入裂纹和设置裂纹面接触，并递交有限元软件Abaqus求解。最后，Franc3D将读取Abaqus的结果文件，以计算裂纹尖端的应力强度因子。
在Abaqus中建立的不含裂纹的平板有限元模型如图3所示，图中不同颜色 区域具有相同的材料属性，但材料名称不同，以便在导入Franc3D时区分裂纹重划分区域。
图3 不含裂纹的平板有限元模型（断裂分析软件Franc3D）
待定义完成后，将上述模型导入到断裂分析软件Franc3D中。需要注意的是，一般情况下，Franc3D会将含有边界条件的单元面或节点设定为保留区域。此时，在植入裂纹或者重划分裂纹的过程中，保留区域的节点布置不会被改变，能够更加方便且精确地将不含裂纹的有限元模型的边界条件映射到含有裂纹的有限元模型的边界条件。如果没有将含有边界条件的单元面或节点设定为保留区域，则这些区域在网格重划分之后可能发生改变，此时Franc3D需要根据原有模型的边界条件，通过外插来重新计算重划分之后模型的边界条件。需要注意的是，无论是否设定保留区域，Franc3D都会执行边界条件的映射操作，但通过设定保留区域，可以尽可能保证映射之后的边界条件仍然是准确的。
但是，在本例中，裂纹重划分区域恰好位于边界条件的施加区域内，在植入裂纹的过程中，必然会破坏原有的节点分布。因此，不能将图3中黄色 区域的上表面设定为保留区域，图4给出了在Franc3D中取消选择保留的边界条件表面的设置，图中Load_Surf即为滚动接触载荷的施加表面。
图4 取消选择保留的边界条件表面
然后，在Franc3D中通过预制的裂纹模版植入穿透型斜裂纹，Franc3D将自动完成裂纹网格的重划分工作。待裂纹网格重划分完成后，在分析菜单中选择执行静态裂纹分析，并勾选定义裂纹面接触，如图5所示。为了与采用Abaqus分析时的模型一致，取接触类型同样为自接触，并且摩擦系数为0.1。
图5 裂纹面接触定义
待定义完成后，即可递交Abaqus执行计算。需要注意的是，与常规计算不同是，本文在施加移动载荷时调用了用户子程序DLOAD，因此在执行计算之前必须指定定义有DLOAD子程序的Fortran文件。一种指定方法为在图6所示的界面中选择查看/编辑命令，并在调用Abaqus的命令行中增加对子程序文件的引用。
另一种方法为在图6所示的界面中勾选写入文件但不运行分析，此时Franc3D将会写入含斜裂纹平板的求解文件（工作文件名+_full.inp），但不会递交Abaqus计算。此时用户需要手动递交Abaqus进行计算。需要注意的是，当使用该方法进行计算时，Franc3D无法直接计算应力强度因子，而是需要用户在计算完成后，在Abaqus后处理界面中运行与求解文件同时生成的writeDtpFile.py脚本文件。该脚本文件将读取Abaqus的数据文件（.odb格式），以提取裂尖区域的位移等参量，并输出为dtp文件。随后，通过在Franc3D中读取该文件，即可计算应力强度因子。
图6 生成求解文件
第二种方法的优点为由于用户可以手动提交Franc3D生成的inp文件，因此对于Franc3D不支持的关键字类型，可以在提交计算之前修改该inp文件，加入Franc3D不支持的关键字类型，然后再提交计算。基于这种方法，理论上可以实现任意复杂载荷条件下的应力强度因子计算。因此，本文采用第二种方法来计算裂纹尖端的应力强度因子，具体的计算流程将在下一节中介绍。
应力强度因子提取方法
应力强度因子提取方法
基于有限元软件Abaqus的提取方法
当采用有限元软件Abaqus进行求解时，可以使用Abaqus内置的围线积分法求解并提取裂纹前缘的应力强度因子。此时，需要在相互作用模块中通过指定裂纹前缘节点和裂纹扩展方向以定义裂纹，并在分析步模块中以时间历程变量的形式定义应力强度因子输出，如图7所示。
图7 应力强度因子输出定义
待计算完成后，裂纹前缘的应力强度因子和节点坐标值将以时间历程变量的形式给出。关于在有限元软件Abaqus中提取应力强度因子的流程，这里不再赘述。
基于断裂分析软件Franc3D的提取方法
当采用Franc3D进行计算时，需要首先提交Franc3D生成的含斜裂纹平板模型的求解文件。需要注意的是，Franc3D会生成三个inp文件，其中以_GLOBAL和_LOCAL结尾的inp文件分别代表重划分区域以外和重划分区域的网格文件，只有以_full结尾的inp文件才是完整的求解文件。
将该求解文件提交Abaqus计算，以生成结果文件（.odb格式）。然后，在Abaqus后处理界面中运行随inp文件同时生成的writeDtpFile.py脚本文件。在本例中，其对应的writeDtpFile.py脚本文件内容如下所示。
from odbAccess import *
import os
# analysis specific data
directory = '.'
file_name = 'Inclined_Crack_CT_full'
node_set = 'LOCAL_CRACKED_NODES'
elem_set = 'TEMPLATE_ELEMENTS'
do_all_frames = True
do_nodal_strain = False
do_integration_sed = False
do_integration_stress = False
# function to output the data for one nodal field value
defOutputNodeValues(fw,dtp_label,frame,odb_key,node_set):
print >> fw, dtp_label
field_output = frame.fieldOutputs[odb_key]
if node_set != None:
field_output = field_output.getSubset(region=node_set)
if len(field_output.values) == 0:
return
data_item = field_output.values[0].data
if field_output.values[0].type == SCALAR:
for field_value in field_output.values:
data = float(field_value.data)
if data != 0.0:
nid = field_value.nodeLabel
print >> fw, nid, float(data)
elif len(data_item) == 3:
for field_value in field_output.values:
nid = field_value.nodeLabel
data = field_value.data
print >> fw,nid,float(data[0]),float(data[1]),float(data[2])
# function to output the data for one element value at the nodes
defOutputElemAtNodeValues(fw,dtp_label,frame,odb_key,elems,elem_set):
print >> fw, dtp_label
field_output = frame.fieldOutputs[odb_key]
if elem_set != None:
field = field_output.getSubset(position=ELEMENT_NODAL,region=elem_set)
else:
field = field_output.getSubset(position=ELEMENT_NODAL)
nc = 0
el_id = 0
for field_value in field.values:
elem = field_value.elementLabel
if el_id != elem:
el_id = elem
nc = 0
else:
nc += 1
conn = elems[elem]
data = field_value.data
if len(data) >= 6:
print >> fw,elem,conn[nc],float(data[0]),float(data[1]),float(data[2]), \
float(data[3]),float(data[4]),float(data[5])
# function to output the data for one element value at the gauss points
defOutputElemAtGpValues(fw,dtp_label,frame,odb_key,elem_set):
print >> fw, dtp_label
field_output = frame.fieldOutputs[odb_key]
if elem_set != None:
field = field_output.getSubset(position=INTEGRATION_POINT,region=elem_set)
else:
field = field_output.getSubset(position=INTEGRATION_POINT)
el_id = 0
ip = 0
for field_value in field.values:
elem = field_value.elementLabel
if el_id != elem:
el_id = elem
ip = 0
else:
ip += 1
data = field_value.data
if field_value.type == SCALAR:
print >> fw, elem,ip,float(data)
elif len(data) >= 6:
print >> fw, elem,ip,float(data[0]),float(data[1]),float(data[2]), \
float(data[3]),float(data[4]),float(data[5])
# function to find element connectivity
defGetElemConn(assembly):
elems = {}
for name, instance in assembly.instances.items():
num_elem = len(instance.elements)
for element in instance.elements:
elems[element.label] = element.connectivity
return elems
# function to check for a uniform temperature
defSameTemp(frame,odb_key):
temp = None
for field_out in frame.fieldOutputs[odb_key].values:
if temp == None:
temp = field_out.data
else:
if field_out.data != temp: returnFalse
returnTrue
# function to process one frame
defDoFrame(frame):
# create a dictionary that maps the first word of a field output
# to the full keys
key_map = {}
for key in frame.fieldOutputs.keys():
v = key.split()
key_map[v[0]] = key
# print >> fw, 'STEPTIME', frame.description
print >> fw, 'TIME', frame.frameValue
# output the data for specific types of data
if key_map.has_key('U'):
OutputNodeValues(fw,'DISPLACEMENT',frame,key_map['U'],node_set)
if key_map.has_key('NT11'):
if SameTemp(frame,key_map['NT11']):
fT = frame.fieldOutputs['NT11'].values[0].data
print >> fw, 'TEMPERATURE'
print >> fw, 'ALL ', fT
else:
OutputNodeValues(fw,'TEMPERATURE',frame,key_map['NT11'],node_set)
if key_map.has_key('CPRESS'):
OutputNodeValues(fw,'PRESSURE',frame,key_map['CPRESS'],node_set)
if key_map.has_key('CSHEAR1'):
OutputNodeValues(fw,'SHEAR STRESS 1',frame,key_map['CSHEAR1'],node_set)
if key_map.has_key('CSHEAR2'):
OutputNodeValues(fw,'SHEAR STRESS 2',frame,key_map['CSHEAR2'],node_set)
if do_nodal_strain:
if key_map.has_key('E'):
OutputElemAtNodeValues(fw,'STRAIN ELEMENT NODAL',frame,key_map['E'],elems,elem_set)
if key_map.has_key('EE'):
OutputElemAtNodeValues(fw,'ELASTIC STRAIN ELEMENT NODAL',frame,key_map['EE'],elems,elem_set)
if key_map.has_key('THE'):
OutputElemAtNodeValues(fw,'THERMAL STRAIN ELEMENT NODAL',frame,key_map['THE'],elems,elem_set)
if do_integration_sed:
if key_map.has_key('SENER'):
OutputElemAtGpValues(fw,'ELASTIC_STRAIN_ENERGY_DENSITY ELEMENT INTEGRATION POINTS',
frame,key_map['SENER'],elem_set)
if key_map.has_key('PENER'):
OutputElemAtGpValues(fw,'PLASTIC_STRAIN_ENERGY_DENSITY ELEMENT INTEGRATION POINTS',
frame,key_map['PENER'],elem_set)
if do_integration_stress:
if key_map.has_key('S'):
OutputElemAtGpValues(fw,'STRESS ELEMENT INTEGRATION POINTS',
frame,key_map['S'],elem_set)
defGetFrameLabel(frame):
# extract the frame 'description' and get the increment
desc = odb.steps[step_name].frames[i].description
di = desc.split()
if di[0] == 'Increment':
ilabel = di[1].split(':')
return ilabel[0]
else:
returnNone
# process the file
os.chdir(directory)
odb = openOdb('%s.odb' % (file_name,))
fw = open('%s.dtp' % (file_name,),'w')
if do_nodal_strain or do_integration_sed or do_integration_stress:
elems = GetElemConn(odb.rootAssembly)
if odb.rootAssembly.instances['PART-1-1'].nodeSets.has_key(node_set):
node_set = odb.rootAssembly.instances['PART-1-1'].nodeSets[node_set]
else:
node_set = None
if odb.rootAssembly.instances['PART-1-1'].elementSets.has_key(elem_set):
elem_set = odb.rootAssembly.instances['PART-1-1'].elementSets[elem_set]
else:
elem_set = None
# loop through the load steps
print >> fw, 'NUM STEP', len(odb.steps.keys())
for step_name in odb.steps.keys():
step = odb.steps[step_name]
step_num = step.number
print >> fw, 'LOADSTEP', step_num
# print >> fw, 'TIMEPERIOD', step.timePeriod
num_frames = len(odb.steps[step_name].frames)
if num_frames == 1:
start_frame = 0
end_frame = 0
else:
start_frame = 0
end_frame = num_frames - 1
ifnot do_nodal_strain: start_frame += 1
ifnot do_all_frames:
start_frame = num_frames - 1
end_frame = start_frame
if end_frame > start_frame:
# substep_cnt = 1
for i in xrange(start_frame,end_frame+1):
ilab = GetFrameLabel(odb.steps[step_name].frames[i])
if ilab != None:
print >> fw,'SUBSTEP',ilab
else:
print >> fw,'SUBSTEP',i
# print >> fw,'SUBSTEP',substep_cnt
DoFrame(odb.steps[step_name].frames[i])
# substep_cnt += 1
else:
DoFrame(odb.steps[step_name].frames[num_frames-1])
fw.close()
可以看到，该脚本文件主要是在Abaqus后处理界面中基于Python语言读取Abaqus生成的脚本文件，并将裂纹尖端节点的位移等信息写入到dtp文件。Franc3D可以读取dtp文件中的信息，并基于这些信息，采用M积分法或其他数值方法计算裂纹尖端的应力强度因子等断裂参量。
在该脚本文件中，directory代表数据文件所在的文件路径，file_name代表数据文件名称，node_set为Franc3D自动创建的节点集，存放了整个裂纹重划分区域的节点，elem_set为Franc3D自动创建的单元集，存放了裂纹前缘模版内的单元，如图8中的红色 区域所示。
图8 裂纹前缘模版内的单元
需要注意的是，如果直接运行该脚本文件，则生成的dpt文件中只会包含结果文件中最后一个增量步的信息。如果需要计算每个增量步下应力强度因子，以输出应力强度因子的时程曲线，则需要将逻辑变量do_all_frames的取值更改为True，即输出每一个增量步下的信息。
待提取完成后，在Franc3D中重新打开Franc3D的项目文件（.fdb格式），并点击文件菜单中的读取结果选项，即可读取脚本文件生成的dpt文件，以计算裂纹前缘的应力强度因子，如图9所示。
图9 dpt文件读取
计算结果分析
计算结果分析
不考虑裂纹面接触的情况
通过取消裂纹面间定义的自接触，图10给出了不考虑裂纹面接触时，由Abaqus和Franc3D计算得到的裂纹前缘中部的I型应力强度因子和II型应力强度因子的变化曲线。
图10 应力强度因子变化曲线（不考虑裂纹面接触）
从图中可以看到，由Abaqus计算得到的I型应力强度因子KI和II型应力强度因子KII与Franc3D给出的计算结果完全吻合，这证明在不考虑裂纹面接触的情况下，基于Abaqus能够准确计算裂纹体的应力强度因子，这符合本文的预期。
考虑裂纹面接触的情况
图11给出了当考虑裂纹面接触时，由Abaqus和Franc3D给出的I型和II型应力强度因子的变化曲线。
图11 应力强度因子变化曲线（考虑裂纹面接触）
从图中可以看到，在考虑裂纹面接触的情况下，I型和II型应力强度因子的峰值要明显偏低，这意味着在不考虑裂纹面接触的情况下，将给出偏于危险的预测结果。同时，由Abaqus给出的II型应力强度因子变化曲线与Franc3D给出的结果较为吻合，而I型应力强度因子的预测结果则存在较大偏差。
为了分析产生偏差的原因，图12给出了在考虑裂纹面接触的情况下，由Abaqus计算得到的不同围线上I型应力强度因子的变化曲线。从图中可以看到，在计算初期，不同围线上的I型应力强度因子非常吻合，这是因为采用围线积分法计算应力强度因子时，应力强度因子的计算值并不依赖于积分路径。然而，当滚动载荷靠近裂纹时，不同围线上的I型应力强度因子出现了明显的发散，对比图11可以看到，这正好对应了Abaqus与Franc3D计算结果存在明显差异的区域。需要注意的是，图11是通过将不同围线上的计算结果取平均值获得的。
图12 不同围线上的I型应力强度因子变化曲线（考虑裂纹面接触）
理论上来说，当裂纹面受到压应力作用而产生闭合时，裂纹尖端的I型应力强度因子应该为零。从图12中可以看到，应力强度因子出现不收敛的起始位置恰好出现在裂纹面产生闭合时。因此，本文认为，Abaqus和Franc3D的计算结果存在差异的原因为：Abaqus在采用围线积分法计算应力强度因子时，并未考虑裂纹面间的接触力对应力强度因子的影响。当滚动载荷远离裂纹尖端时，裂纹面不存在接触力，因此Abaqus与Franc3D的计算结果吻合；而随着滚动载荷靠近裂纹尖端，裂纹面将逐渐接触，并产生接触应力，对于不同围线形成的封闭区域，作用在裂纹面上的接触力显然并不相同，Abaqus并未考虑裂纹面间接触力的影响，导致不同围线上的应力强度因子出现了不收敛的情况。
事实上，在采用围线积分法求解应力强度因子时，通过在围线积分中引入考虑裂纹面接触应力的修正项，可以在考虑裂纹面接触的情况下准确计算裂纹尖端的应力强度因子，但似乎在笔者使用的Abaqus 2020版本以及更高的版本，目前仍未解决该问题。相比之下，由于Franc3D能够考虑裂纹面接触对应力强度因子的影响，因此能够给出更为精确的结果。需要注意的是，从图11中可以看到，当裂纹闭合时，由Franc3D给出的结果要略小于零，这可能是由于Franc3D默认的裂纹模版中，裂尖环形单元的层数较少，导致裂纹尖端区域的接触应力存在较大误差，导致应力强度因子的计算结果不准确，如果增加裂尖环形单元的层数，预计计算结果将更接近零。
此外，在图11中，由Abaqus和Franc3D给出的II型应力强度因子曲线非常吻合，这主要是由于在此种载荷条件下，裂纹面间的摩擦剪切应力较小，因此裂纹面接触对II型应力强度因子的影响不明显。当裂纹面间摩擦剪切应力较大时，不同围线上的II型应力强度因子仍然会出现发散的现象。
总结
总结
本文以受滚动接触载荷作用的含斜裂纹平板为例，分别基于有限元软件Abaqus和断裂分析软件Franc3D求解了裂纹尖端的应力强度因子变化曲线。研究发现，当不考虑裂纹面接触时，Abaqus和Franc3D的计算结果吻合；当考虑裂纹面接触时，Abaqus和Franc3D的计算结果不吻合，Abaqus给出的不同围线上的应力强度因子出现了发散现象，这表明在存在裂纹面接触的情况下，现有版本的Abaqus（Abaqus 2020）不能准确计算裂纹尖端的应力强度因子。
