引言

近年来，中国石油、天然气开发及管道输送得到了迅速发展，管道输送是石油、天然气最经济、合理的运输方式^[1-4]。在石油、天然气等的输送过程中，常常在没有任何预兆的情况下会发生管道腐蚀开裂问题，引发重大的安全事故，导致巨大的经济损失^[5-10]。为保证管道的安全，需要评价管线钢的抗应力腐蚀性能，四点加载法是管线钢抗应力腐蚀性能评价试验方法的一种，见图 1所示。

目前有不少国内外学者开展了大量的钢体材料四点弯曲的试验研究^[11-13]和有限元计算的理论研究^[14-18]。Arun等通过Rousselier模型模拟四点弯曲试验进行了裂纹扩展预测分析^[19]。Aurélien等通过耦合标准的3D数值应用预测了在四点弯曲下环氧树脂/铝试样中的裂纹萌生^[20]。鲜宁等对应力腐蚀试验中四点弯曲加载下带焊缝试样的最大拉应力计算公式进行了详细推导^[21]。

《金属在硫化氢环境中抗硫化物应力开裂和应力腐蚀开裂的实验室试验方法》(GB/T 4157—2017)标准的四点弯曲试样为等厚度的直梁试样^[22]，见图 1所示。而石油行业中的管件，如钻柱、套管、油管以及输送管线等，均为圆形结构管件，标准四点弯曲试样不能真实模拟管件曲梁结构的真实工况。为此，本文设计了一种外壁为完整的一段圆弧，内壁为一直线段和两段对称圆弧构成的非标准、非等厚度的四点弯曲试样，成功地解决了管件四点弯曲试样应力腐蚀开裂的真实工况的模拟问题。又因本文设计试样为非标准试样，标准试样公式不再适用，针对该试样结构推导出了理论公式并与有限元相互验证，确定所推导的理论公式正确性。有限元法计算成本较高，必须要有有限元软件平台才能计算，推导的公式只需简单地把数据代入即可计算出试样外壁任意一点的环向应力，因此，本文设计的试样和推导的理论公式对非标准、非等厚度试样的研究和管件四点弯曲试样的研究都具有很高的实用和推广意义。

1 圆形管件四点弯曲试样结构设计

为了真实模拟管件外壁在最大环向拉应力作用下的应力腐蚀开裂，可以采用C型环和四点弯曲试样来实现。其中，小尺寸管件可直接采用C型环进行加载分析，而对于大尺寸管件，尚无方便方法。因此，本文选取了$R/e$>10的大尺寸试样进行研究($R$—曲率半径，mm；$e$—截面形心到截面内侧边缘距离，mm)，设计并加工出了非标准、非等厚度四点弯曲试样，见图 2所示，并用该试样与中国石油集团石油管工程技术研究院成功地完成了大量的大尺寸管件(外径大于300 mm)应力腐蚀开裂的试验研究工作，为管件曲梁四点弯曲的行业标准化提供了试验依据和理论研究基础。如图 2a所示，该试验装置由非等厚度四点弯曲试样和4个刚性圆柱$B$、$Q$、$B^{'}$和$Q^{'}$组成。该试样外壁$BAB^{'}$为圆弧结构，内壁由两段相等的圆弧$MN$和$M^{'}$$N^{'}$以及直线段$QQ^{'}$组成。整个试样的厚度是变化的，$AP$处厚度$t_0$为最大有效厚度，$G$、$G^{'}$点附近的厚度最小。

该试样的结构特征和优点为：外壁为完整的一段圆弧，内壁为一段直线和两段对称的圆弧构成有效的环向应力变化区和有效的加载区，这样的结构设计有利于推导出试样内外壁精确的几何尺寸的数学计算模型以及其外壁环向应力的力学计算模型。相对标准试样或非标准内外壁均为圆弧的试样，该设计试样更有利于准确加载定位、较小的变形下即可得到较大的环向应力；可通过理论公式和有限元计算同时得到最大环向应力点的位置进行精确标点，更有利于定位和观察最大环向应力点，成功解决了管件四点弯曲真实工况的模拟问题。

确定加载载荷与外壁最大环向拉应力的关系有两种方法，即理论公式计算和有限元法。本文将用这两种方法对非标准四点弯曲试样外壁环向拉应力随其位置路径的变化关系进行研究。

2 试样外壁应力计算理论公式推导

基于试样的对称性，图 2a中非等厚度四点弯曲试样取一半进行理论公式推导。其力学模型和示意图如图 3所示，$O$为全局坐标原点，$OA$为$y$轴，也是对称轴，$B$点为固定不动的刚性圆柱，$Q$点为加载载荷$F$的刚性圆柱，$PQM$为水平直线段，$MNC$为圆弧段，$AEGB$为试样外圆弧段。$\theta$起始点为$y$轴轴线，$OE$为任意$\theta$角的线，其长度始终为试样外半径$d/2$，$OI$为试样内壁半径$R_1$，$R_1$随$\theta$角的变化而变化，其厚度$t$=$OE-OI$=$d/2-R_1$，$t$也随$\theta$角的变化而变化。$O^{'}$为$MN$圆弧段的圆心，也是该圆弧段的局部坐标原点。$C$点为直线$OB$与$MN$圆弧线的交点。

图 3中各变量尺寸和关键点坐标如表 1、表 2所示，尺寸参数与四点弯曲标准试样保持近似的尺寸比，为适合需要可改变试样尺寸。根据图 3的几何结构可得其角度参数的计算式如式(1)$\sim$式(4)，再根据表 1、表 2，可得其基本角度参数的数值如表 3所示，由表 3可知$\theta$角度值为0$\sim$15.698°，$\beta$角度值为50.783°$\sim$82.964°。

$ \begin{align} \left\{ \begin{array}{l} {\theta _1}{\rm{ = }}\dfrac{{180}}{π }{\rm{arctg}} {\dfrac{{{x_Q}}}{{{y_Q}}}} \\[9pt] {\theta _2}{\rm{ = }}\dfrac{{180}}{π }{\rm{arctg}} {\dfrac{{{x_M}}}{{{y_M}}}} \\[9pt] {\theta _3}{\rm{ = }}\dfrac{{180}}{π }{\rm{arctg}} {\dfrac{{{x_N}}}{{{y_N}}}} \\[9pt] {\theta _4}{\rm{ = }}\dfrac{{180}}{π }{\rm{arctg}} {\dfrac{{{x_{{B}}}}}{{{y_B}}}} \end{array} \right. \end{align} $

(1)

$ \begin{align} \left\{ \begin{array}{l} {\alpha _1}{\rm{ = }}\dfrac{{180}}{π }{\rm{arctg}} {\dfrac{b}{a}} \\[9pt] {\alpha _2}{\rm{ = }}\dfrac{{180}}{π }{\rm{arccos}} {\dfrac{{\sqrt {{a^2} + {b^2}} }}{{2{r_1}}}} \\[9pt] \alpha {\rm{ = }}{\alpha _2}{\rm{ - }}{\alpha _1} \end{array} \right. \end{align} $

(2)

$ \begin{align} \left\{\begin{array}{l} {\beta _1}{\rm{ = 2}}(90 - {\alpha _2}) + {\beta _2}\\ {\beta _2}{\rm{ = }}\alpha \\ {\beta _3}{\rm{ = }}\dfrac{{180}}{π }{\rm{arctg}} {\dfrac{{{y_C} - {y_0}}}{{{x_C} - {x_0}}}} \end{array} \right. \end{align} $

(3)

$ \begin{align} \theta = \dfrac{{180}}{π }{\rm{arctg}} {\dfrac{{{x_0} + {r_1}\cos \beta }}{{{y_0} + {r_1}\sin\beta }}} \end{align} $

(4)

根据图 3中的几何关系，可推导出$PM$段的厚度$t$与$\theta$的变化关系式为

$ \begin{align} t = \dfrac{d}{2}\left( {1 - \dfrac{{1 - 2\dfrac{{{t_0}}}{d}}}{{\cos \theta }}} \right)，\theta = {0\sim{\theta _2}} \end{align} $

(5)

式(4)为$MNC$段局部坐标$\beta$角与全局坐标$\theta$角的关系。式(4)的计算比较复杂，基于表 3中$\theta$和$\beta$的数值范围，通过式(4)计算出不同$\theta$下的对应数值$\beta$，并进行拟合，可得$\beta$与$\theta$为线性关系式(6)，其相关系数为0.9999，即完全可以用式(6)代替式(4)。

$ \begin{align} \beta = - 5.079\theta + 130.71，\theta = {{\theta _2} \sim {\theta _4}} \end{align} $

(6)

根据图 3中$MNC$段的圆方程，可推导得该范围内任意$\theta$位置试样的厚度$t$与$\beta$的计算公式为

$ \begin{align} t = \dfrac{d}{2} - \sqrt {{{\left( {{x_0} + {r_1}\cos \beta } \right)}^2} + {{\left( {{y_0} + {r_1}\sin\beta } \right)}^2}} \end{align} $

(7)

由式(5)$\sim$式(7)可知，该非等厚度四点弯曲试样的厚度$t$与角度$\theta$的关系为式(8)。

$ \begin{align} t = f\left( \theta \right)，\theta = 0\sim{\theta _4} \end{align} $

(8)

根据图 3中四点弯曲的力学关系，可推导得$P\longrightarrow$$Q$$\longrightarrow B$内的弯矩变化关系。$PQ$为纯弯矩，不受位置的影响，为一恒定值，$QB$内为线性变化，在支点$B$的弯矩为零，其计算关系见式(9)，式(9)中$x$的计算为式(10)。

$ \begin{align} \left\{ \begin{array}{l} M = F {\dfrac{{H - D}}{2}} , {\rm{ }}\theta = 0\sim{\theta _1}{\rm{}}\\[11pt] M = F{\dfrac{H}{2} - x} , {\rm{ }}\theta = {\theta _1}\sim{\theta _4}{\rm{}} \end{array} \right. \end{align} $

(9)

$ \begin{align} x = \dfrac{d}{2}\sin \theta，\theta={{\theta _1}\sim{\theta _4}} \end{align} $

(10)

曲梁形状的特征可通过曲梁轴线曲率半径$R$与截面形心到截面内侧边缘距离$e$的比值来表示，本文研究的四点弯曲为矩形截面，则$R$刚好等于$(d-t)/2$，$e$等于其厚度的一半，$e=t/2$。对于小曲率梁，即$R/e>10$时，弯曲正应力$\sigma$可以近似用直梁公式(11)计算^[23]，且中性轴通过截面形心，本研究中，$R/e\thickapprox63>10$。

$ \begin{align} \sigma = \dfrac{{My{'}}}{I} \end{align} $

(11)

对四点弯曲矩形截面，中性轴通过形心截面有$y^{'}=t/2$，截面惯性矩$I=wt^3/12$，再由式(9)和式(11)可推导得其截面上外壁的环向拉应力$\sigma$为式(12)。

$ \begin{align} \left\{\begin{array}{l} \sigma {\rm{ = }}\dfrac{{3F(H - D)}}{{w{t^2}}}, {\theta = 0\sim{\theta _1}}\\[11pt] \sigma {\rm{ = }}\dfrac{{3F(H - 2x)}}{{w{t^2}}}, {\theta = {\theta _1}\sim{\theta _4}} \end{array} \right. \end{align} $

(12)

式(12)即为本文推导出的非等厚度四点弯曲外壁环向拉应力$\sigma$随角度$\theta$或试样厚度$t$的变化关系。由于整个试样的厚度$t$是变化的，从图 3中可直观看到，其厚度由大变小，然后由小变大，中间有个最小值厚度$t_{\rm{min}}$。由式(12)可知该最小值$t_{\rm{min}}$必然存在一个最大的环向拉应力$\sigma_{\rm{max}}$。根据式(5)、式(7)和式(9)可分别得出该四点弯曲试样厚度$t$和弯矩$M$随路径角度$\theta$的变化关系，式(9)中弯矩$M$与加载载荷$F$的比值$M/F$仅与试样的结构尺寸有关。由式(5)$\sim$式(10)可得出试样厚度$t$和弯矩$M/F$随路径角度$\theta$的变化关系如图 4和图 5所示，从图 4可见厚度$t$沿路径$\theta$为非线性变化，图 5中弯矩$M/F$为线性变化。

同样，从式(12)可知环向拉应力$\sigma$与加载载荷$F$的比值$\sigma$/$F$仅与试样的结构尺寸有关，并随角度路径$\theta$的变化而变化，计算结果见图 6所示，其$\sigma$/$F$的峰值0.657 8发生在$\theta$=9.373°的位置，只要已知加载载荷的数值大小，该数值乘以图 6中曲线对应的数值，即可得到试样外壁的环向拉应力$\sigma$随角度路径$\theta$的变化关系，或直接用式(12)计算出$\theta$对应的环向拉应力$\sigma$。

3 有限元法与理论公式结果对比分析 3.1 有限元模型建立

前面推导的理论公式(12)只适合于计算试样弹性范围内环向应力计算，不适用于试样发生塑性变形的计算，同时不能计算出试样内部任一点的应力和位移。为了解决这一问题，采用弹塑性有限元方法进行计算，有限元法可得到弹塑性工况下试样内任一点的应力及其应力场分布，为试样提供更完整的应力场数据。

根据本文设计的四点弯曲试样结构可知，属于平面应力带厚度的力学问题，该试样的结构尺寸见表 1。根据图 2中试样的对称性，可建立其有限元力学模型，如图 7所示。在试样的$B$处刚性圆柱施加$x$、$y$方向位移约束，在对称面$AP$处$x$方向位移固定约束，$Q$点刚性圆柱上$x$方向位移约束，在$B$点和$Q$点处的刚性圆柱与试样之间建立接触单元。在$Q$点刚性圆柱上施加一个$y$方向的作用$F$，将在试样外壁上产生环向拉应力$\sigma$，其内壁将产生环向压应力作用，其环向拉应力的最大值发生在环外壁。本文选取石油行业常用管件X56管线钢进行研究分析，其材料的弹性模量为2.1$\times$10$^5$ MPa，泊松比0.3，材料的屈服强度$Y_{\rm{S}}$为386 MPa。

加载过程：首先，在图 7中$Q$刚性圆柱上施加一个很小的载荷$F$，然后，逐渐增加$F$的数值，在该加载过程中，试样外壁的环向拉应力也逐渐增大，直到最大环向拉应力为308.8 MPa，这个数据为其屈服强度$Y_{\rm{S}}$的80%，即可计算出对应的加载载荷$F$，该载荷$F$与$B$处刚性圆柱上的支座反力$F_y$是相等的。如果继续加载到最大环向拉应力为该材料屈服应力386 MPa时，即可计算出试样屈服时的临界加载载荷。

3.2 有限元法与理论公式结果对比

对于X56管线钢某一应力腐蚀开裂试验要求试样最大环向拉应力为80%$Y_{\rm{S}}$，即要求其最大环向拉应力为308.8 MPa。本文建立的四点弯曲的有限元模型如图 7所示，在$Q$点的刚性圆柱上逐渐施加载荷$F$，在试样外壁$AGB$上的环向拉应力也逐渐增大，直到其$G$点附近的最大拉应力为308.8 MPa时，此时对应的载荷$F$值为最终确定的载荷大小。

图 8为有限元计算的试样柱坐标系下环向应力等值线云图，从图 8中可知加载载荷$F$=469.7 N时，试样外壁$G$点的最大环向拉应力$\sigma_{\rm{max}}$=308.8 MPa，即刚好为80%$Y_{\rm{S}}$。从图 8中可见有一明显的中性层，中性层上的应力为零，中性层以上为拉应力，中性层以下为压应力，最大拉应力区发生在$G$点附近区域。

为了详细分析和研究试样外壁$AGB$路径上的环向拉应力的变化规律，将$AGB$圆弧路径上的环向拉应力数据提取出来，其结果见图 9所示，图 9中实线为有限元计算结果曲线，虚线为理论公式(12)的计算结果曲线。从图 9中的曲线结果可知，有限元法的计算结果和理论公式的计算结果非常吻合，即用不同的方法得到了同样的结果，相互验证了其正确性和准确性。理论公式(12)为本文设计的非等厚度的四点弯曲试样装置的标准化提供了可靠的理论计算公式。

为了进一步验证理论公式与有限元法的正确性、可靠性及理论公式的推广应用性，表 4中给出了环向最大拉应力为50%$Y_{\rm{S}}$$\sim$90%$Y_{\rm{S}}$时，5种工况的有限元法和理论公式的计算结果。从表 4中可知，在50%$Y_{\rm{S}}$$\sim$80%$Y_{\rm{S}}$内，理论公式与有限元法的结果误差为$-$0.06%$\sim$$-$0.33%，误差非常小，误差为负值，表明有限元计算结果比理论结果稍偏小。当最大环向拉应力为90%$Y_{\rm{S}}$时，其误差为正的0.13%，表明有限元计算结果比理论结果稍偏大，从有限元法计算结果云图可知，90%$Y_{\rm{S}}$环向拉应力时，其试样内壁$M$点附近开始出现了塑性应力。

综合分析表 4的结果可知，其理论公式(12)的结果与有限元法的结果误差的绝对值为0.06%$\sim$0.33%，几乎没有误差，充分证明了理论公式(12)的正确性、准确性和可推广性。

3.3 外壁$\sigma$/$F$最大值与试样$AP$处厚度$t_0$的关系

由于试样外壁$\sigma$/$F$的最大值(见图 6中峰值)与所研究的最大环向应力直接相关，而试样设计中最重要的参数是试样$AP$处厚度$t_0$，即式(5)中$\theta$=0位置厚度，其位置见图 2所示。在分析中不断改变$t_0$的尺寸数据，由式(12)可以获得每个$t_0$所对应的$\sigma$/$F$最大值，在改变$t_0$的同时，必须保证试样内的最小厚度$t_{\rm{min}}$$\geqslant$3 mm，其结果见图 10所示。根据图 10中的结果曲线可得到最大值$\sigma$/$F$与$t_0$的拟合计算公式见式(13)，其相关系数为0.999 8。

$ \begin{align} \begin{array}{l} {\left. {{\sigma / F}} \right|_{\max }} = 0.1425{t_0}^2 - 2.1799{t_0} + 8.7441{\rm{ }}，\\\hspace{3em} {\rm{ }}{t_0}{\rm{ = 5}}{\rm{.8}} \sim {\rm{7}}{\rm{.0 mm}} \end{array} \end{align} $

(13)

在其他结构尺寸不变的情况下，只要确定式(13)中的试样$AP$处厚度$t_0$，即可得到该结构的$\sigma$/$F$最大值，只要确定其最大环向拉应力$\sigma$的值，即可得到试样试验中的加载载荷$F$的数值。从图 10中曲线的变化规律可知，厚度较小时，即5.8 mm左右其变化斜率较大，且$\sigma$/$F$的数值也较大，因此，要在较小的加载载荷$F$下获得较大的环向应力，厚度越小越好。对于本文研究的非等厚度四点弯曲试样，结合实际应用和试样外壁$\sigma$/$F$最大值随厚度$t_0$的变化关系曲线，建议$t_0$=5.8$\sim$6.5 mm。

4 结论

（1）四点弯曲试样标准为等厚度直梁试样，不能真实模拟管件(如钻柱、套管、油管以及输送管线等)等曲梁结构的真实工况，本文设计了非标准、非等厚度的四点弯曲试样，是一种管件结构四点弯曲应力腐蚀开裂试验的有效方法。

（2）根据非等厚度的四点弯曲试样结构及其受力关系，推导出了其外壁任意一点环向应力与试样加载载荷之间的理论计算公式，为非等厚度的四点弯曲的应力腐蚀开裂试验加载参数提供了简便的计算方法。

（3）文中推导出的理论公式的最大环向应力计算结果与有限元法的计算结果误差的绝对值为0.06%$\sim$0.33%，几乎没有误差，充分证明了理论公式的正确性、准确性和可推广性。

（4）本文设计的非等厚度的四点弯曲试样成功地完成了大尺寸管线应力腐蚀开裂的试验研究工作，为管件曲梁四点弯曲试样的行业试用标准提供了试验依据和理论研究基础，逐渐推荐行业标准或国家标准。

符号说明

$\theta_1$— $OQ$与$y$轴的夹角，(°)；

$\theta_2$— $OM$与$y$轴的夹角，(°)；

$\theta_3$— $ON$与$y$轴的夹角，(°)；

$\theta_4$— $OB$与$y$轴的夹角，(°)；

$x_{\rm{Q}}$— $Q$点$x$坐标，mm；

$x_{\rm{M}}$— $M$点$x$坐标，mm；

$x_{\rm{N}}$— $N$点$x$坐标，mm；

$x_{\rm{B}}$— $B$点$x$坐标，mm；

$x_{\rm{C}}$— $C$点$x$坐标，mm；

$y_{\rm{Q}}$— $Q$点$y$坐标，mm；

$y_{\rm{M}}$— $M$点$y$坐标，mm；

$y_{\rm{N}}$— $N$点$y$坐标，mm；

$y_{\rm{B}}$— $B$点$y$坐标，mm；

$y_{\rm{C}}$— $C$点$y$坐标，mm；

$\alpha_1$— $MN$与$x$轴的夹角，(°)；

$\alpha_2$— $MN$和$O'N$的夹角，(°)；

$\alpha$— $O'N$与$x$轴的夹角，(°)；

$a$— $M、N$两点的$x$坐标差数值，mm；

$b$— $M、N$两点的$y$坐标差数值，mm；

$\beta_1$—局部坐标系中$O'M$与$x$轴的夹角，(°)；

$\beta_2$— $O'N$与$x$轴的夹角，(°)；

$\beta_3$— $O'C$与$x$轴的夹角，(°)；

$x_0$，$y_0$— $O'$点的全局坐标数值，mm；

$r_1$— $MNC$圆弧线的半径，mm；

$\theta$— $OE$与$y$轴的夹角，(°)；

$\beta$— $MNC$段局部坐标角度，(°)；

$t$—试样厚度，mm；

$t_0$—试样$AP$处厚度，mm；

$d$—试样外径，mm；

$H$—外支点间的距离，mm；

$D$—内支点间的距离，mm；

$F$—加载载荷，N；

$M$—弯矩，N$\cdot$mm；

$x$— $QB$内$E$点$x$坐标数值，mm；

$y^{'}$—应力计算点到中性轴的距离，mm；

$I$—截面惯性矩，mm$^4$；

$w$—试样宽度，mm；

$\sigma$—试样环向拉应力，MPa；

$\sigma$/$F$—环向拉应力与加载载荷的比值，mm$^{-2}$。