“貌离神合”的海工结构疲劳分析中的S-N曲线和断裂力学方法
作者 | Simon 仿真秀专栏作者
导读：应用S-N曲线方法分析海洋工程结构物的疲劳问题是目前最为常规的计算手段和设计依据。然而工程实践表明，疲劳破坏案例占到所有结构破坏案例的大多数，远多于屈服和屈曲，这从侧面表明S-N曲线方法可能存在一定缺陷。尽管如此，S-N曲线方法因其直观且易于工程应用的特点，相信今后一段时期内仍然是海洋工程结构物主流的计算分析方法。
我们也可看到近年来，断裂力学方法不断发展。笔者对断裂力学方法在工程上的应用十分关注，目前的主要应用有：
在规范层面，目前船舶行业已经对LNG Type B货舱要求做裂纹扩展分析； 在海工结构（导管架平台）工程应用层面，工程临界分析（ECA）也经常得以应用来分析“已知”裂纹，以支持维修决策和制定检验方案等等, DNV-ST-0119中，对于浮式风机基础，视断裂力学方法为疲劳寿命计算的方法之一。
目前对于疲劳分析方法，应用S-N曲线和断裂力学方法进行分析，无论从书本、规范还是应用都似乎分得很开，有“不相往来”的感觉。
笔者认为研究学习，理解好两者存在的关联，认识断裂力学分析的一些思路和方法对更好得应用S-N曲线方法、一定程度克服其不足很有帮助。本文从工程的角度总结了一些心得体会，抛砖引玉，仅供大家参考。
本文的思路是从大家熟悉的S-N曲线方法入手，讨论应力范围Δσ的意义并引入应力强度因子，建立其与断裂力学方法的联系。再通过一个例子，互验断裂力学方法和S-N曲线方法的结果（附Python代码参考）。主要参考规范DNV-RP-C203以及BS7910。理清一些概念：
应力集中系数SCF和应力强度因子K
ECA和FAD
两种计算方法的联系
对断裂力学，笔者自认为还有很多东西没有掌握，文中难免有误欢迎大家指正。
在采用S-N曲线方法考虑下方左图的普通对接焊细节，一般认为它属于D Curve至F Curve的细节，对应Fatigue Limit (10^7个常应力循环)的应力范围约为40~50 MPa，即应力幅值仅为20至25MPa。给人的直观印象是这个应力水平很低，结构完全处于弹性状态。
需要注意的是S-N曲线中的应力范围Δσ并非“真实”的应力而是“名义”上的应力（或者叫工程应力）。即使采用我们常用的“热点应力”（t x t有限元网格下的应力插值），它也不是真实的应力，而是一种和D曲线“配合”使用的“计算”应力。这可以从相对最接近真实应力的缺口应力（notch stress）来理解。
规范中（DNV-RP-C203）对缺口应力的表达由下图所示。采用非常细密的网格来表达焊接材料和母材的连接区域（半径为1mm的弧线），这样得到的应力需要结合缺口应力S-N曲线（Notch Stress S-N Curve）来计算疲劳损伤。
可以想象要表达半径为1mm的弧线，对单元数量和形状的要求巨大。但其相对比较准确地表达了焊接材料的“几何形状”的影响，故可认为其接近真实的应力水平。
下图中，红笔给出了规范提出的在名义载荷为1时，各缺口应力的值。可知道缺口应力和热点应力之间存在2.5倍的关系（即焊接材料的几何影响）。笔者给出简单的FE分析，得到的结果和规范很相近。
这样来看，对于前面的对接焊细节，10^7循环下疲劳极限的名义应力单幅值为25MPa，对应的缺口应力单幅值可以达到25x1.27x2.5=80MPa。（若不采用设计S-N曲线而采用“平均”S-N曲线，名义应力单幅值为35MPa，对应的缺口应力单幅值可以达到35x1.27x2.5=110MPa）
这个100MPa左右的应力似乎仍然在弹性范围内，但是材料在循环载荷作用下的应力应变曲线和材料本身的拉伸曲线（或者理想弹塑性曲线）是有较大差异的。以32钢为例，如下图所示：
，是σremote，构件几何特征，裂纹长度a的函数。
代入式子
是不是很眼熟呢，对了它就是S-N曲线的表达式。所以从数学表达上来说，S-N曲线可以理解为是断裂力学的“积分形式”，从某种意义上是等效的！
S-N曲线存在Fatigue limit点，即当Δσ很小时，构件不出现疲劳破坏。 这从 ΔKth的角度很好理解，因为只有当ΔKth < ΔK时，裂纹才会扩展。而ΔK正比于Δσ。 S-N曲线和裂纹扩展曲线具有相同的斜率，只是差一个正负符号。 裂纹扩展曲线所分成的I区（裂纹萌生）、II区（裂纹扩展）、III区（裂纹破坏）可以评估得到相对更真实的疲劳损伤情况，比S-N曲线的疲劳损伤线性叠加的概念更合理科学。但同时，采用裂纹扩展曲线需要更深入地考虑载荷循环的先后顺序等因素，工作量有相当的增加。 同S-N曲线一样，近年来裂纹扩展曲线的获得也进行了大量试验。在实际应用中“设计S-N曲线”是将“平均S-N曲线”向“下”移动两个标准差来获得97.7%的“存活概率”（考虑试验数据呈正态分布）。类似的，“设计裂纹扩展曲线”也是将“平均曲线”向“上”平移两个标准差来得到。BS9710推荐的针对焊接节点的单斜率设计裂纹扩展曲线为 
。
Python程序1，计算焊接根部root的疲劳循环次数
# -*- coding: utf-8 -*-"""Verify RP-C203 Sec2.8@author: Simon"""from sympy import *from scipy import integrate # input integrate functionimport numpy as np C = 5.21e-13 # mean plus 2 standard deviationsm = 3.0 # Simplified method for single slope curvesig = 30.0 # G curve stress = 29.24 @ 10^7 cycles M = 1.0B = 25.0 # B = tp =25mm XX = [0.2,0.3,0.4, 0.5, 0.6, 0.7, 0.8, 0.9] # X value for points on curve (2ai/tp) YY = [0.2,0.5,0.64,0.75,0.84,0.93,1.01,1.07] # Y value for points on curve (h/tp)N = [] # Record total Cycles/10^7 def I(a): fwm = (sec(0.5*np.pi*(2*a/W)))**0.5 Mkm = lam0 lam1*(2*a/W) lam2*(2*a/W)**2 Y = M*fwm*Mkm I = ((np.pi*a)**(m/2)*(Y**m))**(-1) return I for i in range(0,8): h = B*YY[i] W = B 2*hlam0 = 0.956-0.343*(h/B)lam1 = -1.219 6.210*(h/B)-12.220*(h/B)**2 9.704*(h/B)**3-2.741*(h/B)**4lam2 = 1.954-7.938*(h/B) 13.299*(h/B)**2-9.541*(h/B)**3 2.513*(h/B)**4 II = integrate.quad(I,B*XX[i]/2,W*0.7/2) # a_initial = 2*ai/2, a_critical = 0.7*W/2 n = II[0]/(C*(sig**m))/(10**7) # n= cycle number / 107 N.append(n)print (N)# The ratio of cycles is about 1.2, which shows the fatigue life of root is longer than the toe.
Python程序2，计算焊接趾端toe的疲劳循环次数，用于对初始缺陷的敏感性分析
# -*- coding: utf-8 -*-"""Sensitive Study on a_initial@author: Simon"""# FM for weld toe on load-carrying cruciform joint import numpy as npfrom sympy import *from scipy import integrate C = 5.21e-13m = 3.0sig = 30.0 # G curve stress = 29.24 @ 10^7 cycles M = 1.0B = 25.0h = 0.5*BL = 2*h B g = B/3 # Typical tapper for P.P.weldingtw = ((B-g)/2 h)*sin(np.pi/4) # P.P. welding throat thickness a_initial = 1 # initial crack size, typically 1mm def Mm(z): return 1.12-0.23*(z/B) 10.6*(z/B)**2-21.7*(z/B)**3 30.4*(z/B)**4 # refer M.4.2 extended surface flaw for a/B <=0.6 def v(z): if z/B <= 0.073: return 0.615*(B/tw)**0.5 else: return 0.83*(B/tw)**0.5 def w(z): if z/B <= 0.073: return -0.31 else: return -0.20def Mk(z): if (v(z)*(z/B)**w(z))>=1.0: return v(z)*(z/B)**w(z) else: return 1.0 def II(a): return (np.pi*a)**(-m/2)*(Mm(a)*Mk(a))**(-1)/(C*sig**m) N, err = integrate.quad(II,a_initial,0.6*B) print (np.log10(N)) # For a_inital = 1mm, total cycle N is about 10^6.93
