受内压航空管道弯曲疲劳可靠性

摘要

引入结合Kriging代理模型法和Monte Carlo抽样法的主动学习可靠性分析方法（AK‒MCS），以管道尺寸、作用载荷以及管材属性为基本随机变量，对受内压悬臂管道在给定工况（即28 MPa内压，自由端5.06 mm纵向位移循环载荷）下的疲劳寿命进行可靠性分析。首先根据管道三点弯曲疲劳试验结果，通过线性回归分析建立管道的Manson‒Coffin疲劳寿命模型，而后采用AK‒MCS方法获得管道的疲劳寿命‒失效概率/可靠度曲线，为内压管道的工程应用提供了基于概率的寿命预测结果。根据该曲线进而得到管道的疲劳寿命频数直方图，发现弯曲疲劳寿命总体上呈正态分布。此外，通过与Monte Carlo方法的比较，证实了AK‒MCS方法在保证受内压航空管道弯曲疲劳可靠性分析精度的同时，还大大降低了计算量。

关键词

内压管道; 疲劳可靠性; Kriging代理模型法; Monte Carlo抽样法; 主动学习方法

航空器中有大量的管道，承担燃油、液压油、润滑油和空气等介质的输送，使飞行器实现如为发动机供油，起落架、襟翼和减速板的收放，润滑降温等工作。管道的失效破裂将会造成严重的飞行事故。随着飞机性能的不断提高，高压力、高功重比、大流量的要求，也导致管道的故障率呈逐年上升态势^［1-2］，因此对管路的失效研究具有非常重要的意义。尽管管道的设计压力是系统工作压力的数倍，如液压管道的设计压力是系统工作压力的6~8倍^［3］，在静载作用下有较大的安全裕度，但在管内流体的高压和外部振动共同作用下会引起疲劳破坏，这是管道故障的主要原因。相较于管内压力，外部振动通过卡箍作用于管道上的动应力会对管道反复作用，并在高应力区引起局部损伤及损伤累积，进而萌生裂纹，当疲劳裂纹长度达到临界值时，裂纹会失稳扩展，导致管道的疲劳破裂失效^［4-7］，因此对受内压管道疲劳寿命的研究是非常必要的。

传统的管道疲劳寿命估计采用的是确定性模型，即在管道寿命分析过程中的参数如管材力学性能、振动载荷及内压、管道尺寸等参数都取确定的数值。但实际上由于加工、装配、测量等原因，上述参数均处于一个不确定的取值范围内，这些参数确定值的选择存在主观性。上述影响因素的随机性导致管道疲劳寿命本质上也是一个随机量。同时，在工程上常采用安全系数法以考虑不确定性因素对管道寿命的影响，但人为给定的安全系数若过大，会导致疲劳寿命的安全裕度过大^［8］，因此，引入了基于概率的可靠性分析方法^［9-10］进行管道疲劳寿命的预测。与确定性分析方法相比，基于概率的可靠性分析方法考虑了上述影响因素的随机性，从概率角度评估了管道寿命，从而可以最大程度上利用管道性能。

有效合理的可靠性分析方法是航空飞行器机结构疲劳寿命设计和控制的基础，一次二阶矩法^［11］（first order second moment method， FOSM）、二次二阶矩法^［12］（second order second moment method， SOSM）、响应面法^［13-14］（response surface method， RSM）以及Monte Carlo法^［15-17］（MCM）等可靠性分析方法已被研究应用。FOSM和SOSM方法的主要思想是在均值点或设计点（失效域内联合概率密度最大的点）处采用线性超平面或者二次超曲面来近似实际功能函数，FOSM方法对于非线性极限状态方程的计算误差较大；利用二次曲面的SOSM可以提高非线性状态方程的计算精度，但由于SOSM需要求解功能函数对随机变量的二阶导数，因此计算很复杂，效率不高^［18-19］。RSM法可通过拟合有限个样本点形成响应面函数，但当求解问题样本点数量大、状态函数高度非线性时，该方法显示出严重精度损失，并且，样本点的选择会影响拟合精度^［20］。MCM具有扎实的理论基础且应用简单，但对于小概率事件，MCM需要进行很多次抽样，导致计算效率较低^［21-22］。Kriging代理模型是目前研究较多的方法^［23-27］，它作为一种估计方差最小的无偏估计模型，具有全局近似与局部随机误差相结合的优点，对非线性程度较高和局部响应突变问题具有良好的拟合效果，是适合用于本问题研究的方法。

本文引入基于概率的可靠性分析方法对受内压悬臂航空管道的弯曲疲劳寿命进行可靠性分析。通过管道三点弯曲试验结果，建立管道的Manson‒Coffin疲劳寿命模型。以管道尺寸、材料参数和载荷为基本随机变量，使用结合Kriging代理模型法和Monte Carlo抽样法的主动学习可靠性分析方法（AK‒MCS）^［28-29］，在管道给定28 MPa内压、自由端5.06 mm纵向位移循环载荷工况条件下，进行管道弯曲疲劳寿命可靠性研究。本文通过与传统MCS方法计算效率的比较，证实了AK‒MCS法可以在保证可靠性分析精度的同时大大降低计算量。另外，根据不同寿命阈值下的失效概率/可靠度的结果，得出管道的疲劳寿命‒失效概率/可靠度曲线，为管道的工程应用提供了基于概率的寿命预测结果。

1 AK‒MCS可靠性分析方法

AK‒MCS可靠性分析方法^［28-29］是Kriging代理模型法和基于Monte Carlo抽样法的主动学习可靠性分析方法，结合了Monte Carlo抽样法和Kriging代理模型法的共同优点，将传统Kriging代理模型的拟合对象从整个空间${\mathbf{R}}^{\mathit{n}}$缩减到仅含有$n_{\mathrm{M}\mathrm{C}}$个样本的Monte Carlo样本池${\mathit{S}}_{\mathrm{M}\mathrm{C}}$上，并利用学习函数在样本池${\mathit{S}}_{\mathrm{M}\mathrm{C}}$中选择对Kriging代理模型拟合贡献最大的下一个评估点更新代理模型，直至建立Kriging代理模型对于样本池${\mathit{S}}_{\mathrm{M}\mathrm{C}}$有着较好的拟合，即自适应学习过程。以下对该方法作简单介绍。

1.1　Kriging代理模型

Kriging代理模型是一种估计方差最小的无偏估计模型，具有局部估计的特点。Kriging代理模型通过样本点（训练点）对模型进行训练以确定随机参数，然后利用建立的Kriging代理模型进行未知响应的预测。Kriging代理模型由参数线性回归模型和高斯随机过程组成：

$$G_{\mathrm{K}}\left(\mathit{x}\right)={\mathit{f}}^{\mathrm{T}}\left(\mathit{x}\right)\mathit{\beta }+z\left(\mathit{x}\right)$$

（1）

式中：$\mathit{f}\left(\mathit{x}\right)={\left[f_{\mathrm{1}}\left(\mathit{x}\right),f_{\mathrm{2}}\left(\mathit{x}\right),\dots ,f_m\left(\mathit{x}\right)\right]}^{\mathrm{T}}$为回归基函数；$\mathit{\beta }={\left[{\beta }_{\mathrm{1}},{\beta }_{\mathrm{2}},\dots ,{\beta }_m\right]}^{\mathrm{T}}$为回归系数；m为基函数的个数；$z\left(\mathit{x}\right)$为一个期望为0，方差为${\sigma }^{\mathrm{2}}$的高斯随机过程，其协方差矩阵定义为

$$\begin{array}{l}\mathrm{C}\mathrm{o}\mathrm{v}\left[z\left({\mathit{x}}_i\right),z\left({\mathit{x}}_j\right)\right]={\sigma }^{\mathrm{2}}R\left({\mathit{x}}_i,{\mathit{x}}_j,\mathit{\theta }\right)\\ i,j=\mathrm{1,2},\dots ,n\end{array}$$

（2）

式中：$R\left({\mathit{x}}_i,{\mathit{x}}_j,\mathit{\theta }\right)$为两个任意的样本${\mathit{x}}_i,{\mathit{x}}_j$之间的相关函数；n为训练样本的个数。相关函数$R\left({\mathit{x}}_i,{\mathit{x}}_j,\mathit{\theta }\right)$有多种形式，常用的是高斯型相关函数，其表达式为

$$R\left({\mathit{x}}_i,{\mathit{x}}_j,\mathit{\theta }\right)=\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \sum _{k=\mathrm{0}}^n}{\theta }_k{\left|x_{ik}-x_{jk}\right|}^{\mathrm{2}}\right)$$

（3）

式中：$x_{ik}$和$x_{jk}$分别是样本向量${\mathit{x}}_i$和${\mathit{x}}_j$的第$k$个分量；${\theta }_k\left(k=\mathrm{1,2},\dots ,n\right)$为相关参数，可以通过极大似然函数取最大值得到最优相关参数${\mathit{\theta }}^{\mathbf{*}}$^［30］，即

$${\mathit{\theta }}^{\mathrm{*}}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}\left\{-\frac{n\mathrm{l}\mathrm{n}{\widehat{\sigma }}^{\mathrm{2}}+\mathrm{l}\mathrm{n}\left|\mathit{R}\right|}{\mathrm{2}}\right\}$$

（4）

式中：过程方差估计值可以通过下式估计，即

$${\widehat{\sigma }}^{\mathrm{2}}=\frac{\mathrm{1}}{n}{\left(\mathit{G}-\mathit{f}\widehat{\mathit{\beta }}\right)}^{\mathrm{T}}{\mathit{R}}^{-\mathrm{1}}\left(\mathit{G}-\mathit{f}\widehat{\mathit{\beta }}\right)$$

（5）

式中：G为训练样本的响应值向量；R为已知训练点的相关矩阵，$R_{ij}=R\left({\mathit{x}}_i,{\mathit{x}}_j,\mathit{\theta }\right),(i,j=\mathrm{1,2},\dots ,n$）；$\widehat{\mathit{\beta }}$为$\mathit{\beta }$的估计值，可以如下式估计，即

$$\widehat{\mathit{\beta }}={\left({\mathit{f}}^{\mathrm{T}}{\mathit{R}}^{-\mathrm{1}}\mathit{f}\right)}^{-\mathrm{1}}{\mathit{f}}^{\mathrm{T}}{\mathit{R}}^{-\mathrm{1}}\mathit{G}$$

（6）

通过求解最优相关参数${\mathit{\theta }}^{\mathbf{*}}$可以得到拟合精度最高的Kriging代理模型。对于一个未知样本点x的响应估计值$G_{\mathrm{K}}\left(\mathit{x}\right)$服从高斯分布，即$G_{\mathrm{K}}\left(\mathit{x}\right)~N\left({\mu }_{G_{\mathrm{K}}}\left(\mathit{x}\right),{\sigma }_{G_{\mathrm{K}}}^{\mathrm{2}}\left(\mathit{x}\right)\right)$，其均值和方差计算式为

$${\mu }_{G_{\mathrm{K}}}\left(\mathit{x}\right)={\mathit{f}}^{\mathrm{T}}\left(\mathit{x}\right)\widehat{\mathit{\beta }}+{\mathit{r}}^{\mathrm{T}}\left(\mathit{x}\right){\mathit{R}}^{-\mathrm{1}}\left(\mathit{G}-\mathit{F}\widehat{\mathit{\beta }}\right)$$

（7）

$${\sigma }_{G_{\mathrm{K}}}^{\mathrm{2}}\left(\mathit{x}\right)={\sigma }^{\mathrm{2}}\left\{\mathrm{1}-{\mathit{r}}^{\mathrm{T}}\left(\mathit{x}\right){\mathit{R}}^{-\mathrm{1}}\mathit{r}\left(\mathit{x}\right)+{\left[{\mathit{f}}^{\mathrm{T}}{\mathit{R}}^{-\mathrm{1}}\mathit{r}\left(\mathit{x}\right)-\mathit{f}\left(\mathit{x}\right)\right]}^{\mathrm{T}}{\left({\mathit{f}}^{\mathrm{T}}{\mathit{R}}^{-\mathrm{1}}\mathit{f}\right)}^{-\mathrm{1}}\left[{\mathit{f}}^{\mathrm{T}}{\mathit{R}}^{-\mathrm{1}}\mathit{r}\left(\mathit{x}\right)-\mathit{f}\left(\mathit{x}\right)\right]\right\}$$

（8）

式中：F为单位列向量；$\mathit{r}\left(\mathit{x}\right)$为训练点和预测点之间的相关函数，$r_i=R\left(\mathit{x},{\mathit{x}}_i,\mathit{\theta }\right),(i=\mathrm{1,2},\dots ,n$）；${\mu }_{G_{\mathrm{K}}}$和${\sigma }_{G_{\mathrm{K}}}^{\mathrm{2}}$可以使用MATLAB的Dace^［31］工具箱进行求解，一般取预测均值${\mu }_{G_{\mathrm{K}}}\left(\mathit{x}\right)$作为Kriging响应$G_{\mathrm{K}}\left(\mathit{x}\right)$，而根据方差${\sigma }_{G_{\mathrm{K}}}^{\mathrm{2}}$的大小，可以判断Kriging代理模型的估计准确程度，这提供了更新Kriging代理模型的指标。

1.2　主动学习可靠性分析方法

结合Kriging代理模型的基于Monte Carlo抽样法的主动学习可靠性分析方法（AK‒MCS）通过以下步骤实现：

（1）Monte Carlo样本池的生成以及初始训练的定义。在输入变量空间${\mathbf{R}}^{\mathit{n}}$根据随机变量的联合概率密度函数生成由$n_{\mathrm{M}\mathrm{C}}$个样本点构成的Monte Carlo样本池并命名为${\mathit{S}}_{\mathrm{M}\mathrm{C}}$。从样本池${\mathit{S}}_{\mathrm{M}\mathrm{C}}$中随机选择$N_{\mathrm{1}}$个随机样本作为初始训练样本，并计算其相应的功能函数响应值G。

（2）Kriging代理模型的建立、预测以及失效概率估计。由选择的$N_{\mathrm{1}}$个样本及相应的响应值建立初始Kriging代理模型，并根据式（1）计算样本池${\mathit{S}}_{\mathrm{M}\mathrm{C}}$所有样本点的功能函数预测值${\mathit{G}}_{\mathrm{K}}$。这一部分主要是在MATLAB的Dace工具箱中实现。然后根据${\mathit{G}}_{\mathrm{K}}$的正负号得到失效概率的估计值${\widehat{P}}_{\mathrm{f}}$，即

$${\widehat{P}}_{\mathrm{f}}=\frac{n_{G_{\mathrm{K}}\le \mathrm{0}}}{n_{\mathrm{M}\mathrm{C}}}$$

（9）

（3）确定性能函数计算的下一个评估点${\mathit{x}}^{\mathit{*}}$并判断学习停止条件。这一阶段需要计算样本池${\mathit{S}}_{\mathrm{M}\mathrm{C}}$中所有样本点的学习函数值，并由此确定下一个评估点${\mathit{x}}^{\mathit{*}}$。现有常用的学习函数有EFF函数^［32］、U函数等^［33］，本文中选择U函数进行下一个评估点${\mathit{x}}^{\mathit{*}}$的选择，其形式如下：

$$U\left(\mathit{x}\right)=\left|\frac{{\mu }_{G_{\mathrm{K}}}\left(\mathit{x}\right)}{{\sigma }_{G_{\mathrm{K}}}\left(\mathit{x}\right)}\right|$$

（10）

$${\mathit{x}}^{\mathit{*}}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}{\mathrm{n}}_{x\in {\mathit{S}}_{\mathrm{M}\mathrm{C}}}U\left(\mathit{x}\right)$$

（11）

（4）通常，$\mathrm{m}\mathrm{i}\mathrm{n}U\left(\mathit{x}\right)\ge \mathrm{2}$可以作为基于U学习函数的Kriging代理模型自适应更新过程收敛的终止条件。当该学习停止条件不满足时，使用下一个评估点${\mathit{x}}^{\mathit{*}}$更新Kriging代理模型。计算下一评估点${\mathit{x}}^{\mathit{*}}$的功能函数值，将这个样本点及其函数值加入之前的设计试验作为新的设计试验，并返回第（2）步重新构建Kriging代理模型。当学习停止条件满足时，继续执行第（5）步。

（5）失效概率变异系数的计算及Monte Carlo样本池的更新。根据式（12）计算失效概率估计值的变异系数，即

$$\mathrm{C}\mathrm{o}{\mathrm{v}}_{\widehat{P_{\mathrm{f}}}}=\sqrt[]{\frac{\mathrm{1}-\widehat{P_{\mathrm{f}}}}{\widehat{P_{\mathrm{f}}}{n}_{\mathrm{M}\mathrm{C}}}}$$

（12）

在自适应学习过程收敛后，通过失效概率的变异系数的大小判断当前失效概率的估计是否稳健。一般变异系数低于5%可以认为样本池${\mathit{S}}_{\mathrm{M}\mathrm{C}}$已经足够大且能够稳健地估计失效概率。若变异系数大于5%，则增加Monte Carlo样本池的数量以提高估计值的稳健性，即扩充样本池${\mathit{S}}_{\mathrm{M}\mathrm{C}}$并返回第（2）步，使用之前建立的Kriging代理模型重新计算得到所有预测值${\mathit{G}}_{\mathrm{K}}^{\mathit{\text{'}}}$，并进行自适应学习，直至满足收敛条件。

（6）结束AK‒MCS方法。当第（5）步计算得到的失效概率估计值的变异系数小于5%时，输出当前计算的失效概率及失效概率估计值的变异系数，完成可靠性评估。

2 管道弯曲疲劳寿命预测

由于航空飞行器中的管道所受绝大多数为弯曲振动，因此本文采用三点弯曲试验方法来确定管道的疲劳寿命。

2.1　管道的三点弯曲疲劳试验

试验参考金属弯曲规范和弯曲疲劳规范^［34-37］，考虑实际卡箍的约束形式，设计夹具如图1所示，实现对称循环疲劳加载，即应力比R = -1。

图1  三点弯曲疲劳试验夹具示意图

Fig.1  Sketch of fixture for three‐point bending fatigue test

在管道的三点弯曲疲劳加载试验中，管道外径12 mm，壁厚1 mm，长度360 mm 。疲劳试验机型号为MTS‒858/2.5T，加载范围±25 kN，加载行程 100 mm，加载频率范围 0~25 Hz。试验选取 5 种位移幅值：2.33、2.67、3.00、3.33和4.00 mm。每种位移幅值保证 3 个以上的有效试验数据。试验在室温下进行，加载频率为18 Hz，数据采样频率为加载频率的10倍。图2为试验照片。

图2  三点弯曲疲劳试验照片

Fig.2  Photo of three‑point bending fatigue test

2.2　疲劳试验数据统计分析

疲劳寿命存在分散性大的问题，主要源自不同样本之间的差异，如材料属性、表面粗糙度和内部缺陷等，因此，需要引入统计理论分析获得应变Δ$\epsilon$与疲劳寿命$N_{\mathrm{f}}$的Δ$\epsilon ‒N_{\mathrm{f}}$曲线。

疲劳试验数据的Δ$\epsilon ‒N_{\mathrm{f}}$关系通过线性回归模型进行统计分析^［38］。该方法假设疲劳寿命服从对数正态分布，并且对数寿命的方差在整个测试范围内是恒定的。在双对数坐标中，Manson‒Coffin疲劳寿命模型可以由$\mathrm{l}\mathrm{g}\left(\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{e}}}{\mathrm{2}}\right)/\mathrm{l}\mathrm{g}\left(\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{p}}}{\mathrm{2}}\right)\mathrm{和}\mathrm{l}\mathrm{g}\left(\mathrm{2}{N}_{\mathrm{f}}\right)$之间的线性模型进行描述，其中$\mathrm{\Delta }{\epsilon }_{\mathrm{e}}$为弹性应变，$\mathrm{\Delta }{\epsilon }_{\mathrm{p}}$为塑性应变，线性关系如图3和图4所示，图中同时给出了线性模型具有95%置信度的置信带。由此可以得到基于疲劳试验数据线性回归分析的Manson‒Coffin疲劳寿命模型为

$$\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{t}}}{\mathrm{2}}=\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{e}}}{\mathrm{2}}$$$$+$$$$\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{p}}}{\mathrm{2}}={\mathrm{10}}^{-\frac{\mathrm{19.540}\mathrm{0}}{\mathrm{8.732}\mathrm{9}}}{\left(\mathrm{2}{N}_{\mathrm{f}}\right)}^{\frac{-\mathrm{1}}{\mathrm{8.732}\mathrm{9}}}+{\mathrm{10}}^{-\frac{\mathrm{2.139}\mathrm{0}}{\mathrm{2.301}\mathrm{3}}}{\left(\mathrm{2}{N}_{\mathrm{f}}\right)}^{\frac{-\mathrm{1}}{\mathrm{2.301}\mathrm{3}}}$$

（13）

图3  $\mathrm{l}\mathrm{g}{N}_{\mathrm{f}}$与$\mathrm{l}\mathrm{g}\left(\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{e}}}{\mathrm{2}}\right)$之间的拟合关系

Fig.3  Fitting relation between $\mathrm{l}\mathrm{g}{N}_{\mathrm{f}}$ and $\mathrm{l}\mathrm{g}\left(\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{e}}}{\mathrm{2}}\right)$

图4  $\mathrm{l}\mathrm{g}{N}_{\mathrm{f}}$与$\mathrm{l}\mathrm{g}\left(\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{p}}}{\mathrm{2}}\right)$之间的拟合关系

Fig.4  Fitting relation between $\mathrm{l}\mathrm{g}{N}_{\mathrm{f}}$ and $\mathrm{l}\mathrm{g}\left(\frac{\mathrm{\Delta }{\epsilon }_{\mathrm{p}}}{\mathrm{2}}\right)$

3 管道弯曲疲劳可靠性分析

以悬臂航空管道在内压28 MPa，自由端5.06 mm纵向位移循环载荷的工况为例^［39］，进行弯曲疲劳寿命可靠性分析。受内压悬臂管道弯曲疲劳试验示意图如图5所示^［39］。该试验结构由直管一端的螺母与阀块之间的螺纹连接，构成悬臂梁的固定端。油泵产生的油液由软管通入阀体为悬臂直管提供固定的内部压强。管道另一端通过堵头密封，防止油液外泄，且在该端通过卡箍将管道与激振器连接。通过激振器的高频上下振动，对管道施加对称的循环位移载荷。

图5  受内压悬臂管道弯曲疲劳试验示意图

Fig.5  Sketch of bending fatigue test of cantilever pipeline under internal pressure

3.1　随机变量的确定

基本随机变量是结构可靠度分析的重要输入变量，其不确定性导致了结构输出响应的不确定性。在受内压管道弯曲疲劳可靠性分析中，基本变量有：几何尺寸、材料参数和工作载荷。这些变量都具有随机性，增加了基本变量选取的复杂性。本文选取的随机变量为：管道的几何尺寸考虑管道半径r和管道厚度t；工作载荷包括位移载荷w和管道内压p；材料参数包括管材的弹性模量$E_{\mathrm{p}}$和泊松比${\mu }_{\mathrm{p}}$以及连接件的弹性模量$E_{\mathrm{c}}$和泊松比${\mu }_{\mathrm{c}}$，以上随机变量及相应的分布模型如表1所示。

表1  基本随机变量及相应的分布模型

Tab.1  Basic random variables and corresponding distribution models

随机变量	单位	均值	方差	分布类型
r	mm	5.000	0.050 0	正态分布
t	mm	0.975	0.009 8	正态分布
w	mm	5.060	0.050 6	正态分布
p	MPa	$28.0$0	0.280 0	正态分布
$E_{\mathrm{p}}$	MPa	$195000$	1950	正态分布
${\mu }_{\mathrm{p}}$		0.280	0.002 8	正态分布
$E_{\mathrm{c}}$	MPa	$200000$	2 000	正态分布
${\mu }_{\mathrm{c}}$		0.300	0.003 0	正态分布

3.2　结构应力分析

本文采用ABAQUS有限元软件进行确定性应力应变分析。根据管道的对称性，为了简化计算，取其1/2的几何模型建立含13 808个单元，19 677个节点的有限元模型，如图6所示。模型采用8节点六面体C3D8R单元进行网格划分，在靠近套筒位置对直管模型进行网格加密。整个航空管道内壁施加28 MPa内压，另外在加载端施加幅值为5.06 mm的循环弯曲位移载荷。计算悬臂管道的弯曲应力应变，并由式（1）计算相应的疲劳寿命。在确定性分析中，各参数均选择表1所示的随机变量的均值进行分析计算。

图6  受内压管道有限元模型

Fig.6  FEM model of internal pressured pipeline

确定性分析的节点最大主应变和主应力结果如图7和图8所示。由图7和图8可知，管道的最大主应力和主应变出现在固定端一侧的管道根部，即此位置的疲劳寿命最低。由图7可得，此位置的半应变幅$\mathrm{\Delta }{\epsilon }_{\mathrm{t}}/\mathrm{2}=\mathrm{3.290}\mathrm{8}\times {\mathrm{10}}^{-\mathrm{3}}$，由公式（1）计算得到$N_{\mathrm{f}}=\mathrm{12}\mathrm{031}$个周期，该寿命与Yang等^［39］试验测得的受内压管道疲劳平均寿命$N_{\mathrm{f}}=\mathrm{12}\mathrm{806}$非常接近，验证了该寿命模型的准确性和有效性。同时，记录单次有限元计算时间为47 s，用以后续的AK‒MCS法和MCS法的计算效率对比。

图7  悬臂管主应变分布云图和最大主应变位置局部放大图

Fig.7  Nodal principal strain distribution contour of cantilever pipeline and local amplified view of maximum principal strain position

图8  悬臂管主应力分布云图和最大主应力位置局部放大图

Fig.8  Nodal principal stress distribution contour of cantilever pipeline and local amplified view of maximum principal stress position

3.3　受内压管道疲劳可靠性分析及结果

根据确定性分析结果，受内压管道的最大应力应变位置出现在固定端一侧的管道根部，该位置的疲劳寿命即为管道的寿命。本文以受内压悬臂管道在固定端一侧的管道根部位置的疲劳寿命小于等于规定的寿命阈值$N_{\mathrm{f}}^{\mathrm{*}}$作为失效事件F，以该部位的疲劳寿命大于$N_{\mathrm{f}}^{\mathrm{*}}$作为安全事件S。由此可以定义受内压管道疲劳寿命可靠性分析的极限状态函数为

$$G\left(\mathit{X}\right)=N_{\mathrm{f}}\left(g\left(\mathit{X}\right)\right)-N_{\mathrm{f}}^{\mathrm{*}}$$

（14）

式中：$N_{\mathrm{f}}\left(·\right)$为疲劳寿命函数；$g\left(·\right)$为结构响应计算函数，描述输入随机变量与输出的应力应变之间的关系；$N_{\mathrm{f}}^{\mathrm{*}}$为寿命阈值，不同的寿命阈值对应不同的失效概率$P_{\mathrm{f}}$。相应的管道疲劳的失效概率和相应的可靠概率定义如下：

$$P_{\mathrm{f}}=P\left\{F\right\}=P\left\{G\left(\mathit{X}\right)\le N_{\mathrm{f}}^{*}\right\}$$

（15）

$$P_{\mathrm{s}}=P\left\{S\right\}=P\left\{G\left(\mathit{X}\right)>N_{\mathrm{f}}^{*}\right\}$$

（16）

式中：P｛｝表示概率算子。

本文采用AK‒MCS分析方法根据2.2节中的步骤对受内压管道弯曲疲劳的可靠性进行分析。在步骤（1）中，由表1中随机变量的联合概率密度函数生成$n_{\mathrm{M}\mathrm{C}}$=50 000个Monte Carlo样本点，并从中随机选择$N_{\mathrm{1}}$=30个随机样本作为初始训练样本s（s₁，s₂，…，s₃₀）。在步骤（2）中，通过ABAQUS软件计算得到训练样本对应的σ_m和 Δε_t ，根据式（13）计算相应的疲劳寿命N_f ，结合式（14）计算寿命阈值$N_{\mathrm{f}}^{\mathrm{*}}$对应的功能函数G值。

为了验证AK‒MCS法的有效性，本文同时采用Monte Carlo法进行结果比对。由于Monte Carlo法在求解小概率事件时存在所需样本点过大的问题，故针对疲劳寿命阈值为9 500个周期的小失效概率事件进行效率对比。两种方法的效率对比结果见表2。对于Monte Carlo法，进行50 000次随机抽样并计算其功能函数的值。根据4.2节的结构应力分析可知，每次有限元计算时间为47 s，所以共需要花费652.78 h，计算量很大。而相应使用代理模型的方法，共需要调用141次有限元计算，其中包括30次建立初始Kriging模型和111次修正模型，共花费1.84 h，降低了99.72%的有限元计算时间，而调用Kriging代理模型计算样本点响应的时间相较于有限元计算时间可以忽略不计。并且两种方法得到的失效概率基本上可以认为是相同的。故AK‒MCS方法对于失效概率求解有着很好的准确性和高效性。

表2  MCS法和AK‒MCS法计算效率对比

Tab.2  Comparison of computation efficiency between MCS and AK‒MCS

可靠性分析方法	失效概率$/\mathrm{\%}$	有限元调用次数	计算时间/h
MCS	0. 081	50 000	652.78
AK‒MCS	0. 082	30+111=141	1.84

取不同的寿命阈值$N_{\mathrm{f}}^{\mathrm{*}}$，可以得到相应的失效概率和可靠度，结果如表3所示。

表3  受内压管道疲劳可靠性分析结果

Tab.3  Reliability analysis results of fatigue life of pressure pipeline

$N_{\mathrm{f}}^{\mathrm{*}}$	9 500	9 600	9 700	9 800	9 900	10 000	10 100	10 200	10 300	10 400	10 500	10 600	10 700
$P_{\mathrm{f}}/\mathrm{\%}$	0.08	0.13	0.20	0.31	0.52	0.78	1.14	1.60	2.20	3.12	4.20	5.34	6.83
$P_{\mathrm{s}}/\mathrm{\%}$	99.92	99.87	99.80	99.69	99.48	99.22	98.86	98.40	97.80	96.88	95.80	94.66	93.17
$N_{\mathrm{f}}^{\mathrm{*}}$	10 800	10 900	11 000	11 100	11 200	11 300	11 400	11 500	11 600	11 700	11 800	11 900	12 000
$P_{\mathrm{f}}/\mathrm{\%}$	8.53	10.54	12.86	15.28	18.00	20.81	23.83	27.29	31.00	35.03	39.37	43.64	48.02
$P_{\mathrm{s}}/\mathrm{\%}$	91.47	89.46	87.14	84.72	82.00	79.19	76.17	72.71	69.00	64.97	60.63	56.36	51.98
$N_{\mathrm{f}}^{\mathrm{*}}$	12 031	12 100	12 200	12 300	12 400	12 500	12 600	12 700	12 800	12 900	13 000	13 100	13 200
$P_{\mathrm{f}}/\mathrm{\%}$	49.89	52.18	56.32	60.43	64.77	68.60	71.96	75.24	78.09	80.55	82.87	85.42	87.60
$P_{\mathrm{s}}/\mathrm{\%}$	50.11	47.82	43.68	39.57	35.23	31.40	28.04	24.76	21.91	19.45	17.13	14.58	12.40
$N_{\mathrm{f}}^{\mathrm{*}}$	13 300	13 400	13 500	13 600	13 700	13 800	13 900	14 000	14 100	14 200	14 300	14 400	14 500
$P_{\mathrm{f}}/\mathrm{\%}$	89.50	91.33	92.74	94.01	95.20	96.26	97.35	98.28	98.57	98.83	99.00	99.16	99.26
$P_{\mathrm{s}}/\mathrm{\%}$	10.50	8.67	7.26	5.99	4.80	3.74	2.65	1.72	1.43	1.17	1.00	0.84	0.74

疲劳寿命‒失效概率/可靠度曲线如图9所示，该曲线可以为内压管道的工程应用提供基于概率的寿命预测结果。

图9  疲劳寿命‒失效概率/可靠度曲线

Fig.9  Fatigue life-failure probability/reliability curve

由表3和图9可知，受内压管道的疲劳寿命小于10 000个周期的概率在1%以内，而当疲劳寿命达到14 500个周期时，受内压管道失效的概率将超过99%。同时，确定性分析寿命$N_{\mathrm{f}}$对应的可靠度为49.89%。

根据图9中的疲劳寿命‒失效概率曲线，可以得到样本总数为50 000的情况下，受内压管道疲劳寿命频数直方图，如图10所示。由图10曲线可知，受内压管道疲劳寿命分布总体上为正态分布形式，曲线峰值出现于疲劳寿命均值12 035处。受内压管道的疲劳寿命分布可以通过一定量的疲劳试验样本均值和方差进行近似估计，该结论对受内压管道的初步工程设计具有指导意义。

图10  受内压管道疲劳寿命频数直方图

Fig.10  Histogram of fatigue life frequency of pipeline under internal pressure

4 结论

在高流体内压和外部随机振动共同作用下的疲劳失效是飞行器管道的主要破坏模式之一，因此相应的疲劳寿命的研究是非常必要的。本文首先针对航空管道进行无内压三点弯曲疲劳试验，试验结果表明，管道的疲劳寿命存在着很大的分散性，这种分散性来自于不同样本之间的差异，如材料属性、表面粗糙度和内部缺陷等。故通过线性回归模型进行统计分析，建立基于统计理论的受内压管道Manson‒Coffin疲劳寿命模型；其次，研究受内压管道在给定工况（28 MPa内压，自由端5.06 mm纵向位移循环载荷）的疲劳寿命，以管道尺寸、作用载荷以及管材材料性能参数为基本随机变量；最后，采用AK‒MCS方法进行受内压管道疲劳寿命可靠性分析。分析得到以下结论：

（1）由AK‒MCS方法得到受内压管道的弯曲疲劳寿命‒失效概率/可靠度曲线，该曲线为受内压管道的工程应用提供了基于概率的寿命预测值。根据该曲线得到样本总数为50 000下的受内压管道疲劳寿命的频数直方图，疲劳寿命总体上呈正态分布。

（2）在管道弯曲疲劳寿命可靠性分析过程中，通过确定性分析计算得到各变量取均值情况下的疲劳寿命为$N_{\mathrm{f}}=\mathrm{12}\mathrm{031}$，该寿命与YANG等人试验测得的受内压管道疲劳平均寿命$N_{\mathrm{f}}=\mathrm{12}\mathrm{806}$非常接近，从而验证了Manson‒Coffin疲劳寿命模型的有效性和准确性。同时，确定性分析寿命$N_{\mathrm{f}}$对应的可靠性为49.89%，从而验证了AK‒MCS法可靠性分析方法的有效性和准确性。

（3）通过AK‒MCS方法与传统的MCS方法计算失效效率的对比发现，对于寿命阈值取9 500个周期的小概率失效事件，AK‒MCS方法将调用有限元计算次数从传统的MCS方法需要调用的50 000次降低到141次，降低了99.72%的有限元计算时间，同时保证计算所得的失效概率几乎相等。从而证明AK‒MCS可以在保证可靠性分析精度的同时大大降低计算量。

作者贡献声明

沈兴铿：试验实施，理论分析和计算，论文撰写。

王光强：试验方法指导。

杨 婧：试验指导。

员婉莹：可靠性分析理论指导。

杨宏伟：试验设计和实施。

张屹尚：试验设计和指导。

戴 瑛：论文总体指导和修改。

贺鹏飞：方法和试验指导。