隧道在交通运输领域发挥着至关重要的作用，然而衬砌劣化所引发的坍塌等事故对国民经济以及人员生命安全产生了重大影响. 1960年伦敦地铁北线隧道由硫酸盐侵蚀引起了衬砌破损，整个维修工作持续了6 a并耗费1 500万英镑^[1]. 1990年，香港地铁隧道在仅运行20 a后发生了钢筋锈蚀，使用超高压射水技术才完成了隧道衬砌的更新^[2]. 2011年河南的南阳隧道以及2012年重庆的九龙隧道也都由于衬砌劣化发生坍塌，造成了恶劣的影响. 世界各地由隧道衬砌劣化引起经济损失的事例层出不穷^[3-6]. 合理预估衬砌的劣化过程对于隧道养护以及指导衬砌设计具有重要的意义.

衬砌的抗弯承载能力是表征隧道抵御外部荷载、保证安全运营的重要参数之一. 目前对于衬砌的抗弯承载能力已有大量研究. Abbas等^[7]采用全尺寸试验研究衬砌在纯弯曲以及循环荷载下的抗弯承载能力. Liu等^[8-11]通过一系列全尺寸试验研究管片衬砌在不同工况荷载下的极限抗弯承载能力. Liu等^[12-13]采用理论与试验分析的方法给出纤维混凝土的抗弯承载能力计算模型. 这些研究为隧道衬砌极限抗弯承载能力的研究提供了基础，但是其均着眼于衬砌的短期承载能力，忽略了衬砌受环境影响在长期服役中发生的性能劣化.

相关研究表明，由钢筋锈蚀所引起的开裂、保护层脱落、钢筋截面积损失以及黏结力下降对衬砌的承载能力有极大的影响^[14]. 在材料层次，Saetta等^[15-16]结合理论与模拟，预测钢筋混凝土由于碳化侵蚀发生钢筋锈蚀的初始时间；Almusallam^[17]采用试验方法研究钢筋锈蚀程度对其力学性能的影响；Coronelli^[18]通过试验对钢筋锈蚀程度对其与混泥土的黏结力的影响进行深入研究. 这些针对钢筋混凝土在材料层次劣化的研究为后续构件层次劣化模型的构建提供了基础. 许多学者针对钢筋截面损失、保护层开裂、黏结力下降对钢筋混凝土梁承载能力的影响进行了深入研究. 然而，由于运营环境不同，利用此类模型描述隧道衬砌的劣化仍存在困难. 除此之外，影响隧道衬砌长期性能劣化的因素还包括植根材料本身性能参数的不确定性、人为施工造成构件尺寸的不确定性以及环境变量随机作用的不确定性等. 大量的工程实践证明这些不确定性会对构件承载能力造成影响^[19-26]，不过目前鲜有针对考虑不确定性的隧道衬砌承载能力的计算模型的研究.

本研究结合有限元衬砌外部荷载计算、时变承载能力概率模型建立以及蒙特卡洛模拟（Monte-Carlo simulation，MCS）方法结果所服从的分布分析构建隧道衬砌可靠度劣化模型的分析框架. 1）通过研究得到时变不确定性的衬砌抗弯承载力模型. 相关试验表明纵向钢筋截面积的减少以及钢筋混凝土之间黏结力的下降是造成衬砌性能劣化的主要原因^{[15-16, 18, 27-30]}. 基于此，建立同时考虑钢筋截面面积减少以及黏结力下降引起的变形不协调的衬砌混凝土双筋截面模型，考虑受压区混凝土达到强度极限建立衬砌时变抗弯承载力计算模型. 2）通过对时变抗弯承载力模型中参变量的分析，结合统计学原理建立考虑参变量不确定性的抗弯承载能力时变概率模型. 3）采用MCS方法验证提出的概率模型的准确性以及计算效率. 通过分析MCS方法结果，对承载力模型输出结果所服从的分布进行研究. 4）通过有限元软件计算衬砌所受荷载，结合得到的时变承载力概率模型构建可靠度分析功能函数，采用结构安全联合委员会推荐（joint committee，JC）方法对衬砌服役期的可靠率与失效概率进行详细研究. 得到衬砌可靠度劣化的规律，以期为隧道维护策略的制定以及衬砌设计提供理论支持.

大量研究给出了经典的单筋截面梁承载能力的计算方法，本研究在此基础之上建立适用于隧道衬砌的双筋截面模型，并通过钢筋锈蚀率反映模型的时变特性.

试验表明钢筋截面减小以及钢筋混凝土黏结退化是钢筋混凝土结构抗弯承载能力减弱的主要因素^[14]. 钢筋截面减小程度以及钢筋-混凝土黏结力下降导致的变形不协调程度均可以用钢筋锈蚀率（锈蚀部分质量占总质量的比重）来表示. 假设钢筋发生均匀锈蚀，锈蚀率的表达式为

(1) $\eta (t) = 1 - {{{{\left( {{d_{\rm{b}}} - {d_{\rm{e}}}} \right)}^2}} / {{d_{\rm{b}}}^2}}.$

式中：η（t）为t时刻钢筋的锈蚀率；d_b为钢筋直径；d_e为锈蚀深度. d_e的表达式^[14]为

(2) ${d_{\rm{e}}} = 0.011\;6Jt.$

其中，J为腐蚀电流密度，t为时间.

考虑到锈蚀引起的钢筋截面减小，受拉区锈蚀钢筋的拉力表达式为

(3) $T\left( t \right) = \left[ {1 - {\eta _{\rm{t}}}\left( t \right)} \right]{A_{\rm{s}}}{\sigma _{\rm{s}}} = \left[ {1 - {\eta _{\rm{t}}}\left( t \right)} \right]{A_{\rm{s}}}{E_{\rm{s}}}{\varepsilon _{\rm{s}}}.$

式中：A_s为受拉区钢筋总截面积，E_s 为钢筋弹性模量，σ_s为钢筋应力，ε_s为受拉钢筋应变，η_t(t)为t时刻受拉区钢筋的截面锈蚀率.

受压区钢筋的合力表达式为

(4) $T'(t) = \left[ {1 - {\eta _{\rm{c}}}\left( t \right)} \right]{A'_{\rm{s}}}{\sigma _{\rm{s}}} = \left[ {1 - {\eta _{\rm{c}}}\left( t \right)} \right]{A_{\rm{s}}}{E_{\rm{s}}}{\varepsilon '_{\rm{s}}}.$

式中： ${A'_{\rm{s}}}$为受压区钢筋总截面积； ${\varepsilon '_{\rm{s}}}$为受压区钢筋应变；η_c(t)为t时刻受压区钢筋的截面锈蚀率.

在锈蚀钢筋混凝土构件中，由于钢筋和混凝土会发生相对滑移，钢筋和混凝土之间不能满足变形协调几何关系. 由文献[31]的试验结果可以发现，截面上的混凝土平均应变仍然符合平截面假定. 因此，截面曲率 $\phi $与混凝土的应变存在如下关系：

(5) $\phi = \frac{{{\varepsilon _{{\rm{cu}}}}}}{{{x_{\rm{c}}}}} = \frac{{{{\varepsilon _{\rm{c}} '}}}}{{{x_{\rm{c}}} - {{a_{{{\rm{s}}}} '}}}}.$

式中：ε_cu为混凝土的极限压应变，x_c为中性轴高度， ${\varepsilon '_{\rm{c}}}$为混凝土压应变， ${a'_{\rm{s}}}$为受压区混凝土边缘到钢筋外缘的距离.

用等效的矩形压应力（见图1（d））代替实际的压应力分布（见图1（c））. 等效矩形由参数β、γ确定. β为等效矩形高度与中和轴高度的比值，γ为等效矩形的应力与受压区应力极值的比值，表达式分别为

(6) $\beta = {x / {{x_{\rm{c}}}}},\; \gamma ={\sigma / {{\sigma _{\rm{0}}}}}.$

式中：x为受压区高度，σ为等效矩形的应力，σ₀为受压区应力极值.

将 $\;\beta = {x / {{x_{\rm{c}}}}}$代入式（5），可以得到

(7) ${\varepsilon '_{\rm{c}}} = {{\beta {\varepsilon _{{\rm{cu}}}}(x/\beta - {{a_{{s}}'}})} / x}.$

锈蚀钢筋混凝土的钢筋和混凝土在同一点处的不协调程度表达式^[31]为

(8) $m\left( {\eta \left( t \right)} \right) = {{{\varepsilon _{\rm{c}}}} / {{\varepsilon _{\rm{s}}}}}.$

式中：ε_c为混凝土应变.

徐善华^[31]通过试验得到 $m(\eta )$的表达式为

(9) $m\left( {\eta \left( t \right)} \right)=0.043\;8\eta {\left( t \right)^2} + 0.051\;3\eta \left( t \right) + 0.993\;6.$

将式（8）代入式（7），可以得到

(10) ${\varepsilon '_{\rm{s}}} = \frac{{\beta {\varepsilon _{{\rm{cu}}}}(x/\beta - {{a_{{s}}'}})}}{{xm\left( {\eta \left( t \right)} \right)}}.$

由图1（b）可知

(11) $ \left. {\begin{aligned} {\frac{{{\varepsilon _{{\rm{cu}}}}}}{{{\varepsilon _{{\rm{cs}}}}}}=\frac{{{x_{\rm{c}}}}}{{{h_0} - {x_{\rm{c}}}}} = \frac{{x/\beta }}{{{h_0} - x/\beta }} = \frac{{{\varepsilon _{{\rm{cu}}}}}}{{m\left( {\eta \left( t \right)} \right){\varepsilon _{\rm{s}}}}}}, \\ {\frac{{{\varepsilon _{{\rm{cu}}}}}}{{{\varepsilon '_{\rm{s}}}}}=\frac{{{x_{\rm{c}}}}}{{{x_{\rm{c}}} - {a_1}}}} .\qquad\qquad\qquad \end{aligned}} \right\} $

式中：ε_cs为锈蚀钢筋的应变，h₀为截面有效高度，a₁为受拉区钢筋中心距保护层外缘距离.

假设在衬砌破坏时受压区混凝土达到极限应变ε_cu. 由 $C + {\rm{ }}T' = T$，可得

(12) ${f_{{\rm{cd}}}}bx + \left[ {1 - {\eta _{\rm{c}}}\left( t \right)} \right]{A'_{\rm{s}}}{E_{\rm{s}}}{\varepsilon '_{\rm{s}}} = \left[ {1 - {\eta _{\rm{t}}}\left( t \right)} \right]{A_{\rm{s}}}{E_{\rm{s}}}{\varepsilon _{\rm{s}}}.$

式中：f_cd为混凝土抗压强度；C为混凝土所受压力合力； $T' $为受压区钢筋压力合力；T为受压区钢筋拉力合力；b为截面宽度，取单位长度.

将式（6）、（11）代入式（12），可以得到

(13) $\begin{split} & {f_{{\rm{cd}}}}bm{\eta _{\rm{t}}}\left( t \right){x^2} + \left\{ {\left[ {1 - {\eta _{\rm{t}}}\left( t \right)} \right]{A_{\rm{s}}}{E_{\rm{s}}}{\varepsilon _{{\rm{cu}}}} + }\right.\\ & \quad \left. {\left[ {1 - {\eta _{\rm{c}}}\left( t \right)} \right]{{A_{\rm{s}}'}}{{E_{\rm{s}}'}}{\varepsilon _{{\rm{cu}}}}m{\eta _{\rm{c}}}\left( t \right)} \right\}x - \\ & \quad \left[ {1 - {\eta _{\rm{c}}}\left( t \right)} \right]{{A_{\rm{s}}'}}{{E_{\rm{s}}'}}m\left( {{\eta _{\rm{t}}}\left( t \right)} \right)\beta {\varepsilon _{{\rm{cu}}}}{a_{\rm{s}}} - \\ & \quad\left[ {1 - {\eta _{\rm{t}}}\left( t \right)} \right]{A_{\rm{s}}}{E_{\rm{s}}}\beta {\varepsilon _{{\rm{cu}}}}{h_0} = 0. \end{split}\qquad\quad $

式中： ${{E_{\rm{s}}'}} $为受压区钢筋弹性模量，a_s为保护层厚度.

考虑到受拉区钢筋黏结力退化，在衬砌破坏时无法确定其是否已到达抗拉强度极限，受拉区钢筋合力点处的弯矩表达式为

(14) $\begin{split} {M_{\rm{u}}} =& {f_{{\rm{cd}}}}bx\left( {{h_0} - x/2} \right) + \\ & {{\left[ {1 - {\eta _{\rm{c}}}\left( t \right)} \right]{{A_{\rm{s}}'}}{{E_{\rm{s}}'}\;}\beta {\;\varepsilon _{{\rm{cu}}}}(x/\beta - {{a_{{s}}'}})({h_0} - a{'_{\rm{s}}})} / x}. \end{split} $

将所提出公式的计算结果与文献[32]结果进行对比，验证公式的准确性. 文献[32]试验的模型如图2所示，模型参数及试验与理论结果的对比情况如表1所示. 表中，f_sd为钢筋屈服强度，η为钢筋锈蚀率，P_m为试验梁破坏时所施加的荷载，P_s为使用本研究模型计算得到的试验梁破坏荷载. 对比发现3个试件计算结果的误差分别为0.02、0.26、0.04，平均误差为0.11，计算公式具有较高的准确性.

参数的不确定性对模型输出结果有较大的影响. 由式（13）、（14）可以看出，衬砌抗弯承载能力与混凝土抗压强度、钢筋锈蚀率等参数相关.

抗弯承载能力的计算模型可以简写为

(15) ${M_{\rm{u}}} = f({{X}},t),$

(16) ${{X}} = \left[{f_{{\rm{cd}}}},{A_{\rm{s}}},A{'_{\rm{s}}},{a_{\rm{s}}},a{'_{\rm{s}}},h,{E_{\rm{s}}},{{E_{\rm{s}}'}},{d_{\rm{b}}},{J_{\rm{c}}},{J_{\rm{t}}} \right] .$

式中：J_c为受压区钢筋腐蚀电流密度；J_t为受拉区钢筋腐蚀电流密度，h为截面高度.

考虑到衬砌一般为对称配筋，式（16）可以简化为

(17) ${{X}} = \left[ {{f_{{\rm{cd}}}},{A_{\rm{s}}},{a_{\rm{s}}},h,{E_{\rm{s}}},{d_{\rm{b}}},{J_{\rm{c}}},{J_{\rm{t}}}} \right].$

根据统计学原理，抗弯承载力概率模型的均值可以表示为

(18) $\begin{split} E({M_{\rm{u}}}) = & {\mu _{{M_{\rm{u}}}}} = f({{{\mu}} _{{X}}},t)= \\ & f({\mu _{{f_{{\rm{cd}}}}}},{\mu _{{A_{\rm{s}}}}},{\mu _{{a_{\rm{s}}}}},{\mu _h},{\mu _{{E_{\rm{s}}}}},{\mu _{{d_{\rm{b}}}}},{\mu _{{J_{\rm{c}}}}},{\mu _{{J_{\rm{t}}}}},t). \\ \end{split} $

方差可以表示为

(19) $ {\rm{Var}}({M_{\rm{u}}}) = {\sigma ^2}_{{M_{\rm{u}}}} = \mathop \sum \limits_{i = 1}^n {\left\{ {{{\left.{\frac{{\partial f({X},t)}}{{\partial {X}}}} \right|}_{{\mu X},t}}} \right\}^2}{\sigma_{{X}}}. $

在得到隧道的初始计算参数以及地层等模型之后，衬砌所受的荷载可以通过有限元（finite element method, FEM）软件计算. 在得到衬砌承载弯矩后，结构可靠度指标计算的功能函数可以表达为

(20) $Z = {g_X}({{X}}) = {M_{\rm{u}}} - M.$

式中：M_u的均值和方差可分别由式（18）、（19）求出，M为外部荷载作用下的承受弯矩.

可靠度指标表达式为

(21) $\beta =\frac{{{g_X}({{{x}}^*}) + \displaystyle\sum\limits_{i = 1}^n {\left( {{{\partial {g_X}({{{x}}^*})} / {\partial {X_i}}}} \right)\left( {{\mu _{{X_i}}} - {x_i}^ * } \right)} }}{{{{\left( {\displaystyle\sum\limits_{i = 1}^n {{{\left[ {{{\partial {g_X}({{{x}}^*})} / {\partial {X_i}}}} \right]}^2}\sigma _{{X_i}}^2} } \right)}^{1/2}}}}.$

式中： ${{{x}}^ * }={\left[ {{x_1}^ * ,{x_2}^ * , \cdot \cdot \cdot ,{x_n}^ * } \right]^{\rm{T}}}$为极限状态面上的点， $\sigma _{{X_i}}^2 $为变量X_i的方差.

定义变量X_i的灵敏度系数为

(22) ${\alpha _{{X_i}}}=\cos\; {\theta _{{X_i}}}=\frac{{\left( {{{\partial {g_X}({{{x}}^*})}/ {\partial {X_i}}}} \right)\sigma _{{X_i}}^{}}}{{{{\left( {\displaystyle\sum\limits_{i = 1}^n {{{\left[ {{{\partial {g_X}({{{x}}^*})} / {\partial {X_i}}}} \right]}^2}\sigma _{{X_i}}^2} } \right)}^{1/2}}}}.$

极限状态面上点的坐标未知，假设各参变量均值点为初始验算点 ${{{x}}^*}=\left( {{\mu _{{X_1}}},{\mu _{{X_2}}}, \cdot \cdot \cdot ,{\mu _{{X_n}}}} \right)$. 验算点在原始空间中的坐标为

(23) ${x_i}^ * = {\mu _{{X_i}}} + \beta {\sigma _{{X_i}}}\cos \;{\theta _{{X_i}}}.$

当存在非正态变量时，须对其进行当量正态化后再带入迭代计算. 设非正态变量X={X₁，X₂， $\cdots $，X_i， $\cdots $，X_n}中元素的均值为 ${\mu _{{X_i}}}$，标准差为 ${\sigma _{{X_i}}}$，概率密度函数为 ${f_{{X_i}}}\left( {{x_i}} \right)$，累积分布函数为 ${F_{{X_i}}}\left( {{x_i}} \right)$. 与X_i对应的当量正态化变量 ${X_i}^\prime $的均值和标准差分别为 ${\mu _{{{X_i'}}}}$、 ${\sigma _{{{X_i'}}}}$，概率密度函数与累积分布函数分别为 ${f_{{{X_i'}}}}\left( {{{x_i'}}} \right) $、 ${F_{{{X_i'}}}}\left( {{{x_i'}}} \right)$. 当量正态化的条件为在验算点 ${{{x}}^ * }$处X_i、 ${X_i}^\prime $的累积分布函数和概率密度函数分别对应相等，即

(24) ${F_{{{X_i'}}}}\left( {{x_i^ *} } \right)=\varPhi \left( {\frac{{{x_i^ *} - \mu _{{X_i'}}}}{{{\sigma _{{{X_i'}}}}}}} \right)={F_{{X_i}}}\left( {{x_i^ *} } \right),$

(25) ${f_{{{X_i'}}}}\left( {{x_i^ *} } \right) = \frac{1}{{{\sigma _{{{X_i'}}}}}}\varphi \left( {\frac{{{x_i^ *} - \mu_{{X_i'}}}}{{{\sigma _{{{X_i'}}}}}}} \right) = {f_{{X_i}}}\left( {{x_i^ *} } \right).$

式中：Φ为标准正态分布的累积分布函数，φ为标准正态分布的概率密度函数。

当量正态化变量的均值与标准差表达式分别为

(26) ${\mu _{{{X_i'}}}}={x_i}^ * - {\varPhi ^{ - 1}}\left[ {{F_{{X_i}}}\left( {{x_i}^ * } \right)} \right]{\sigma _{{{X_i'}}}},$

(27) ${\sigma _{{{X_i'}}}}={{\varphi \left\{ {{\varPhi ^{ - 1}}\left[ {{F_{{X_i}}}\left( {{x_i}^ * } \right)} \right]} \right\}} {\Big /} {{f_{{X_i}}}\left( {{x_i}^ * } \right)}}.$

衬砌劣化可靠度计算的具体计算框架如图3所示. 以功能函数的求解为节点，整个计算的流程可以分为两部分. 在节点之前为功能函数的求解过程，主要包括3个步骤：1）抗弯承载能力时变概率模型的求解，通过式（18）、（19）求解得到任意时刻衬砌抗弯承载力的均值与方差；2）利用MCS 方法获取抗弯承载能力分布；3）利用有限元软件获取衬砌荷载模型. 通过第1）、2）步的求解可以获得可靠度计算的抗力模型，结合第3）步获取的荷载模型可以得到可靠度计算的功能函数，从而进行节点后的JC法可靠度迭代计算过程. 对不同时间节点的可靠度进行循环计算，可以得到隧道衬砌寿命周期内的可靠度劣化情况. 须指出的是，在由迭代过程求解出可靠度指标β_k后，在按式（23）求解新的验算点 ${{{x}}^*}$时，须区分正态变量与非正态变量的 ${\mu _{{X_i}}}$. 当参数服从正态分布时，在每次循环时 ${\mu _{{X_i}}}$取参数的均值；当参数服从非正态变量时，由于须进行当量正态化， ${\mu _{{X_i}}}$可以直接参与循环计算.

基本参数的概率统计特征是衬砌抗弯承载力时变概率特征分析的基础. 结构安全联合委员会（joint committee on structural safety，JCSS）给出了大量钢筋、混凝土等基本参数的概率统计特征值，其余参数的统计特征可通过长期收集或参考其他文献得到. 在本研究中，为了验证所述方法的应用及准确性，钢筋、混凝土以及截面几何参数的变异系数均取JCSS推荐值. 考虑到隧道内侧为大气环境，假定内侧的腐蚀电流密度为0. 以某深埋公路隧道为例，隧道的衬砌配筋如图4所示.

围岩级别为IV级，衬砌重度为25 kN/m³，围岩容重为20 kN/m³，围岩弹性模量为5 GPa，衬砌泊松比为0.25，围岩泊松比为0.30，参照JTGD70-2004《公路隧道设计规范》^[33]，假设衬砌承担围岩荷载的40%. 考虑地下水中氯离子对衬砌的侵蚀，假定衬砌与围岩接触一侧腐蚀电流密度J_out= 0.5 μA/cm^2[4]，内侧不受氯离子侵蚀影响，腐蚀电流密度J_in= 0. 其余计算参数如表2所示. 须说明的是在隧道衬砌中一般同时存在正负弯矩（衬砌内拉外压为正），受压钢筋与受拉钢筋腐蚀电流的表达式分别为

(28) $\left. \begin{array}{l} {{\text{正弯矩}}}:\;\;\;\;{J_{\rm{c}}} = {J_{{\rm{out}}}},\;{J_{\rm{t}}} = {J_{{\rm{in}}}},\\ {{\text{负弯矩}}}:\;\;\;\;{J_{\rm{c}}} = {J_{{\rm{in}}}},\;{J_{\rm{t}}} = {J_{{\rm{out}}}}. \end{array}\right\}$

如图5所示为第2章建立的理论模型与1 000 000次MCS的衬砌抗弯承载力计算结果的对比. 图中，μ、σ分别为均值、标准差. 可以看出，理论模型结果与MCS方法的计算结果较一致，均值计算的最大误差为3.5%，标准差计算的最大误差为5.0%. 电脑配置为CPU CORE i5，Memory 16 GB；理论计算时长为115 s. MCS时长为9 400 s. 在理论计算过程中，内存占用为0.48 GB，MCS方法的内存占用为2.84 GB.

由图5可以看出，在正弯矩区域由于受拉侧钢筋受到腐蚀，抗弯承载能力降低幅度较大. 如图6、7所示，对MCS方法求解的0、100 a的承载力模拟数据进行分布拟合. 图中，P为概率，P_c为累积概率. 由图可知，在初始时刻模拟数据与正态分布以及对数正态分布的拟合结果均较接近；随着时间的增加，正态分布的拟合结果与模拟数据之间的偏差变大，但是对数正态分布拟合结果与模拟数据的分布仍较吻合，只是在概率曲线的两端（模拟数据出现可能性小）存在偏差. 可以认为，随着时间的发展衬砌抗弯承载能力服从对数正态分布，则在任意时刻衬砌抗弯承载能力的概率密度函数与累积分布函数表达式分别为

(29) $f(x,t) = \frac{1}{{{\sigma _{{{M_{\rm{u}}}}(t)}}x\sqrt {2\pi } }}\exp\; \left[ { - \frac{1}{2}{{\left( {\frac{{\ln\;x - {\mu _{{{M_{\rm{u}}}}(t)}}}}{{{\sigma _{{{M_{\rm{u}}}}(t)}}}}} \right)}^2}} \right],$

(30) $F(x,t) = \frac{1}{2} + \frac{1}{2}{\rm{erf}}\;\left[ {\frac{{\ln x - {\mu _{{{M_{\rm{u}}}}(t)}}}}{{{\sigma _{{M_{\rm{u}}}(t)}}\sqrt 2 }}} \right].$

按照式（29）、（30）绘制不同时间节点的抗弯承载能力的概率密度函数（probability density function，PDF）与累积分布函数（cumulative distribution function，CDF），如图8、9所示. 随着时间的发展，抗弯承载力的概率密度函数与累积分布函数向左移动，即抗弯承载力均值随时间减小. 此外，概率密度函数由平坦变集中，表示抗弯承载力方差随时间减小. 从第20 a到第100 a，曲线越来越密集，说明均值与方差在腐蚀后期变化缓慢. 另外两图的累积分布函数的范围以及曲线斜率也能反映这一规律. 对比两图可知，正弯矩区域内的抗弯承载力受时间影响较大.

计算模型取5倍洞径范围，采用反转应力荷载释放法. 该方法使用地层结构法计算衬砌受力，具体操作为：1）计算整个地层开挖隧道前的初始地应力；2）开挖隧道，激活衬砌位置单元的混凝土材料属性，与此同时，对衬砌单元施加反向的60%的初始地应力，以模拟隧道的自承能力；3）利用FEM软件计算得到内力结果（见图10）；4）按照图2流程在MATLAB中编写计算程序，计算得出可靠度与时间、角度的关系三维曲面，如图11所示. 图中，θ为不同位置衬砌到拱顶的角度.

由图10可以看出，衬砌自拱顶至拱肩范围内承受正弯矩，拱肩至拱脚范围内承受负弯矩，仰拱范围内承受正弯矩. 拱顶与仰拱的正弯矩较大，分别为23.5、42.7 kN·m. 最大负弯矩出现在拱脚处，为88.8 kN·m. 由图11中可靠度曲面在θ-o-β平面的投影可以看出，正弯矩作用区域与负弯矩作用区域的可靠指标有较大差异，负弯矩区域可靠度高于正弯矩区域. 负弯矩作用区域可靠指标最小值远高于3.7（《建筑可靠度统一标准》^[34]中一级结构构件的目标可靠度指标参考值为3.7）. 由三维曲面在t-o-β平面的投影可以看出，随着时间的变化衬砌的抗弯承载可靠度指标变化较大，在前20 a可靠度指标降低较为迅速，之后降低较为缓慢.

考虑最不利状况，拱顶向下160°处的可靠度指标与对应的失效概率随时间的变化关系如图12所示. 图中，P_f为失效概率可以看出，在衬砌钢筋发生锈蚀后不到30 a的时间内，其可靠度指标已不足3.7，对应失效概率已高于0.01%，由此可见，推后钢筋开始锈蚀的时间对保障隧道在服役期内满足目标可靠度具有重要意义.

（1）随着时间的增长，抗弯承载力的均值和标准差均减小；

（2）当初始计算变量同时存在正态分布变量与对数正态分布变量时，抗弯承载力计算结果更服从对数正态分布；

（3）衬砌截面不同位置处的可靠度指标不同，在研究中应当注意最危险截面.

（4）研究考虑了隧道受地下水中氯离子侵蚀的抗弯承载特征. 其他腐蚀因素对衬砌可靠性影响也可以通过本研究框架实现.

（5）对弯矩正负不同区段须区别考虑拉压钢筋对应的锈蚀率.

（6）衬砌腐蚀电流密度与隧道环境温度、湿度等因素相关，接下来拟就这些因素展开研究.