索杆结构[1-2]是一类应用广泛的空间结构，而结构中的钢索构件在长期服役过程中内部金属会产生位错滑移、原子扩散和晶界滑动，其宏观现象表征为蠕变和应力松弛现象[3-4](本文统一表述为蠕变)．受其影响，结构中钢索的长度将随时间逐步增加，结构的预应力水平会逐步降低，进而导致结构刚度的削减和性能的退化．文献[5]进行了不同初始应力下小直径高强钢丝的松弛试验并考虑了预应力混凝土构件中收缩徐变影响的变长度松弛试验结果，但在索杆体系中的适用性有待进一步验证．文献[6-7]用核函数的乘积和两种加性近似，给出了钢索在不变荷载下的非线性蠕变一般本构方程，但该方程中的蠕变系数依赖于应力水平，而实验中的应力水平是不连续的，采用不连续应力水平下的蠕变系数对结构进行分析必然会造成误差．文献[8]提出了简单的直线钢绞线的有限元模型，并发现钢丝受轴向荷载的响应是由钢丝间的扭动和滑动控制的，这对于蠕变机理的描述具有重要意义．文献[9]利用蠕变松弛理论和弹性理论，建立了钢丝的应力松弛模型，并对1×7钢绞线在拉伸荷载下的应力松弛进行了研究．文献[10]对高强钢丝进行了高温蠕变试验，并建立了高强钢丝应力松弛的计算模型，然而该模型仅适用于高温情况，无法考虑常温钢丝的蠕变．文献[11]研究了锻造压机机架的预应力钢丝的蠕变规律，由常温试验得到了两种应力下140 h的蠕变数据并提出了预应力损失模型，但是其试验时间较短，无法用于结构钢索的蠕变预测中．以上研究大多针对构件材料的蠕变属性进行分析，对于索杆结构整体性能受蠕变影响的分析则较少．且在实际工程中，结构安全性能的鉴定仍旧是基于结构现状的力学性能进行判断的，通过结构的健康监测及预警功能确保结构性能安全的研究仍在发展中；同时，国内的相关规范也处于探索阶段，结构在剩余服役期间的安全性能难以被准确评估．本研究考虑了索杆结构的蠕变效应，建立了相应的承载能力极限状态与正常使用极限状态下的极限状态方程与不等式．通过蠕变发展模型及力学推导，给出了不等式方程中参数的求解方法．同时基于数学规划并以钢索蠕变随时间的非线性变化规律作为约束条件，提出了确定结构蠕变伸长限值及安全工作年限的具体方法．最后在算例中对两种极限状态下不同可靠度指标对应的蠕变伸长限值及安全工作年限进行了分析，并与规范限值进行了对比，可为结构的可靠性及实际工程中的监测预警提供参考．1 索杆结构的两种极限状态分析在建筑结构的使用过程中，结构是以可靠(安全、适用、耐久)和失效两种状态存在的．两者状态的界限就称之为结构的极限状态，其实质是结构工作状态的一种阈值，若超过这一阈值，则结构处于不安全或者不耐久、不适用的状态；反之则结构处于可靠状态．我国《建筑结构可靠性统一标准》[12]中将结构的状态分为承载能力极限状态与正常使用极限状态．对于索杆结构，蠕变效应会导致钢索构件长度增加及结构刚度的削减，主要表征为钢索发生松弛现象，节点位移增加．当钢索的预应力小于某一限值时，钢索将退出工作，结构无法维持正常的形状．因此，结构的承载能力极限状态应考虑钢索在蠕变作用下的最小应力问题．而节点位移的变化也应符合相应规范要求，此类问题应归属于正常使用极限状态．结构的状态可以用结构状态方程表示为Z(G1,G2,…,Gn)=m0+∑i=1nmiGi，(1)式中：Z为结构功能函数；Gi为独立随机变量(i=1，2，…，n)；m0为影响Z的常量；mi为随机变量的权重系数．结构的安全可靠性能则可以通过可靠度指标β来衡量[13]，β值越大表示结构越安全，反之则表示结构的失效概率越大，β=m0+∑i=1nmiμGi/∑i=1n(miλGi)2，(2)式中μGi和λGi分别为变量Gi的均值和标准差．规范中规定了β在两种极限状态下的取值，然而实际工程中对于结构可靠度指标β的取值往往还须结合工程经验综合判定．这主要是由于可靠度指标取值本身难以通过实验方法确定，规范对于β的确定仍在探索中，因此对其进行理论上的分析是必要的．以下从索杆结构的两种极限状态具体分析．1.1　承载能力极限状态分析处于工作状态下的索杆结构的构件预应力水平较高，结构的性能变化近似线性，因此不同索的蠕变伸长对于同一根索的内力及节点位移的影响可以线性叠加．设结构共有ns个钢索，nn个节点，各个钢索之间发生蠕变伸长是相互独立的．考虑到索杆结构多为对称结构，由蠕变效应所引起的应力下降往往会在对称位置的同类构件中同步发生．当单根钢索的预应力降低至某一限值退出工作时，其对称位置的同类构件也几近失效或已同时失效．且索杆结构体系冗余度较低，某一类构件失效时极易导致结构无法继续维持正常形状，因此可以认为此时结构将突破承载能力极限状态．建立第i根钢索的极限状态方程Zc=Ti-αciTsi，(3)式中：Ti为考虑蠕变效应后的第i根钢索的内力值；Tsi为第i根钢索维持结构外观及刚度的最小预张力；αci为Tsi的安全系数，其取值越大，钢索的最小预张力储备越多，当αci=1时，结构及构件将不再考虑额外的最小预张力储备．又有Ti=T¯i+∑j=1nsmijΔLj，(4)式中：T¯i为第i根钢索不考虑蠕变影响的内力值；ΔLj为第j根钢索的蠕变伸长量；mij为第j根钢索发生单位蠕变伸长时对第i根钢索内力的变化参数．Ts可按下式计算[14]    Ts≥W[(L/4)cos θ+(3/4)wsin θ]/(2w)≈20Wcos θ,(5)式中：W为拉索自重；L为拉索跨度；w为不影响结构外观的拉索垂度，一般取w≤L/160；θ为拉索与地面的夹角．根据式(2)和(4)建立第i根钢索达到极限状态时的可靠度指标βi方程βi=T¯i-αciTsi+∑j=1nsmijμj/∑j=1ns(mijλj)2，(6)式中μj和λj分别为第j根构件蠕变伸长量ΔLj的均值和标准差．为了保证结构的安全可靠性，须使得式(6)中βi大于限值β¯，即β≥β¯．代入式(6)，化简可得到以下极限状态不等式∑j=1nsmijλj2≤T¯i-αciTsi+∑j=1nsmijμj/β¯2．(7)将以上不等式进行整合并写成矩阵形式，可得Mnsλ≤Tsβ，(8)式中：Mns定义为钢索内力变化参数矩阵，Mns=m112m122⋯m1ns2m212m222⋯m2ns2⋮mns12mns22⋯mnsns2；λ为蠕变伸长标准差矩阵，λ=[λ12,λ22,…,λns2]T；Tsβ为钢索内力限值矩阵      Tsβ=T¯1-αc1Ts1+∑j=1nsm1jμj/β¯2,T¯2-αc2Ts2+∑j=1nsm2jμj/β¯2,…,                T¯ns-αcnsTsns+∑j=1nsmnsjμj/β¯2T.1.2　正常使用极限状态分析当节点变形超过某一限值时，结构将突破正常使用极限状态．建立第k个节点的极限状态方程Zd=αdiHck-Hk，(9)式中：Hck为第k个节点的容许变形限值，本研究取结构跨度的1/250进行控制；Hk为考虑蠕变伸长影响后的第k个节点构件的变形值；αdi为Hck的安全系数，其取值越小，容许变形储备越多，当αdi=1时，结构及构件将不再考虑额外的容许变形．又有Hk=H¯k+∑j=1nsηkjΔLj，(10)式中：H¯k为第k个节点不考虑蠕变影响的变形值；ηkj为第j根钢索发生单位蠕变伸长时对第k个节点变形的变化参数．同理建立第k个节点达到极限状态时的可靠度指标βi=αdiHck-H¯k-∑j=1nsηkjμj/∑j=1nsηkjλj2(11)及极限状态不等式∑j=1nsηkjλj2≤αdiHck-H¯k-∑j=1nsηkjμj/β¯2．(12)整理以上极限状态不等式，并将其表示成矩阵形式，可得Nnnλ≤Hcβ，(13)式中：Nnn定义为节点变形变化参数矩阵，Nnn=η112η122⋯η1ns2η212η222⋯η2ns2⋮ηnn12ηnn22⋯ηnnns2；Hcβ为节点变形限值矩阵，    Hcβ=αd1Hc1-H¯1-∑j=1nsη1jμj/β¯2,αd2Hc2-H¯2-∑j=1nsη2jμj/β¯2,…,αdnnHcnn-H¯nn-∑j=1nsηnnjμj/β¯2T.2 极限状态不等式求解由式(2)可知，可靠度指标βi是变量Gi的均值μGi与标准差λGi的多元函数，对其求解前先须确定这两类极限状态不等式的参数矩阵，即式(8)中的Mns和Tsβ与式(13)中的Nnn和Hcβ．2.1　变化参数矩阵确定索杆结构的节点数为nn，钢索数为ns，其静力平衡方程可以表达为Kx=A1q=P,(14)式中：K为刚度矩阵；x为节点坐标矢量；A1为对应于力密度q的平衡矩阵；P为外荷载矢量．基于小变形假定，结构中构件的内力向量为t=(EAL0-1)×V=(EAL0-1)×(L-L0),(15)式中：E为构件弹性模量对角矩阵；A为构件截面面积矢量；L0和L分别为构件初始长度与变形后的长度矢量；V为变形矢量．变形后构件长度对节点位移的偏导等于协调矩阵∂L∂x=∂(L0+V)∂x=∂V∂x=BL,(16)且BL与平衡矩阵A1具有如下关系BLT=A1L-1.(17)构件力密度向量q可以表达为q=L-1t=EAL0-1-L-1.(18)将式(18)代入平衡方程式(15)中可以得到Kx=A1EAL0-1-L-1=P.(19)结构的力学性能受到变量x和L的影响，因此在式(14)中对x和L进行全微分，可以得到(∂Kx/∂x)dx+(∂Kx/∂L)dL=dP,(20)式中∂Kx/∂x=KT，KT为结构的切线刚度矩阵．将式(17)和(19)代入式(20)并整理可得dx=KT-1dP-KT-1BLTCndL，(21)式中Cn=EAL0-1为构件刚度的对角矩阵．式(21)中由钢索蠕变伸长所影响的节点位移变化部分可由dxL=-KT-1BLTCndL来表示．当dL取单位长度时，dxL矩阵中第k行第j列元素的物理意义为第j根构件发生单位蠕变伸长时对第k个节点变形的影响，即式(10)中的变化参数ηkj．对其进行平方并集合为矩阵形式，即可得到节点变形变化参数矩阵Nnn．同理，式(15)对x和L进行全微分，可以得到dt=Cn(BLdx+dL).(22)将式(21)代入式(22)可以得到新的dt表达式dt=CnBLKT-1dP+Cn[Ins-CnBLKT-1BLT]dL.(23)式(23)中由钢索蠕变伸长所影响的内力变化部分可由dtL=Cn[Ins-CnBLKT-1BLT]dL.同理，dtL矩阵中的第i行第j列元素的物理意义为第j根构件发生单位蠕变伸长时对第i根构件内力的影响，即式(4)中的变化参数mij．对其进行平方并集合为矩阵形式，即可得到构件内力变化参数矩阵Mns．2.2　钢索的蠕变发展模型由于极限状态不等式含有两类变量μ和λ，因此无法直接求解变量．为此引入钢索的蠕变发展模型，来确定各个钢索的蠕变均值．文献[15-16]开展了常温下钢索的长期蠕变实验，并基于黏弹性理论模型拟合了蠕变的本构方程，通过Ansys建立符合实际蠕变性能的有限元模型，验证了蠕变本构方程的准确性．本研究采用该模型进行蠕变伸长ΔLi的计算，模型如下    ε1=ε2(100)×t/100    (t≤100 h);    ε2=(-2.04×10-4σ0.113d0.2+4.44×10-4)·exp[(0.007 88σ0.076d-0.013 7-0.013 68)t]+0.001 95× σ0.056d0.086-0.002 9   (100 ht4 300 h);    ε3=(9.097×10-22σ4.752d-1.203+1.809×10-9×(t-4 300)+ε2(4 300)    (t≥4 300 h),(24)式中：d为钢索直径；σ为钢索的应力；t为时间．该模型是通过长期蠕变实验数据拟合而成，根据大数定律可以认为该模型整体误差较小，接近钢索蠕变的实际情况，因此可以用于式(6)和式(11)中的蠕变伸长量均值μj的计算．此时βi仅为λ1,λ2,…,λns的函数．当给定可靠度指标限值时，即可以确定λ1,λ2,…,λns的范围．根据拉依达准则[17]，蠕变伸长限值与其均值和标准差存在以下关系dLi=μi+3λi.(25)根据式(25)可以计算出钢索的蠕变伸长限值，结合蠕变发展模型式(24)即可以对结构的安全年限进行预测计算．可见λ1,λ2,…,λns的求解是计算钢索蠕变伸长限值与结构安全年限的关键．2.3　索长容许限值的数学规划分析以上推导将βi转化为λ1,λ2,…,λns的多元函数，此类函数的求解属于数学规划方法．合理的目标函数与可行的边界约束是求解数学规划问题的基础．将式(8)和(13)作为可靠度边界的非线性约束，其物理意义为可靠度限值控制不等式在设计空间中形成的超曲面，其内部的空间为可行域．每一个可行域中的设计点均满足约束条件，且设计点距离原点越近，可靠度越高，反之则越低．为了计算出满足结构安全前提下的最大蠕变伸长限值，须要在可行域中找出距离原点最远的设计点，因此该问题为单目标非线性规划问题．其设计变量和目标函数如下：λ=[λ1,λ2,…,λi,…,λns]T；(26)min  f1(λ)=min-∑i=1nsλi2,(27)式中f1(λ)为钢索蠕变伸长限值标准差和的负值，该函数越小表明构件蠕变伸长限值越大．式(24)中ε为关于时间t、直径d和钢索应力σ的函数．根据式(24)计算出服役年限te时刻下不同钢索的蠕变伸长ΔLi并将其归一化后，得到ΔL¯1,ΔL¯2,…,ΔL¯ns，将其作为条件约束建立蠕变约束方程，同时设计变量λ的数学意义为一个非负数，结合以上条件建立约束方程如下 λ1ΔL¯1=λ2ΔL¯2=⋯=λnsΔL¯ns    (λi≥0)．(28)式(28)在可行域中为一直线，该直线上的点均满足年限te时刻下蠕变比例，其与约束界面的交点即是满足结构安全符合te时刻下蠕变比例的设计点，根据式(25)即可计算出te年份下的蠕变伸长限值dLi．须要注意的是，蠕变伸长限值dLi与时间t无关，仅由目标函数和约束方程控制，由此计算得出的dLi可能与实际蠕变伸长量ΔLi呈倍数关系，且蠕变应变受结构刚度变化影响，随时间呈非线性变化，不同年限下各索的蠕变比例并不相同，即式(28)会随时间而变化．因此，应根据dLi与ΔLi值的大小调整te，重新计算dLi，直至两者相等．考虑到结构寿命较长，式(24)的蠕变都是以年小时数为单位进行计算的，不同年份下ΔLi的结果是离散的，因此与ΔLi差值绝对值最小的dLi即可认为是实际的蠕变伸长限值．建立约束方程s.t.  min∑insΔLi-dLi,(29)满足约束方程式(29)的dLi其对应的年限te即为安全年限．当服役时间tte时，钢索的蠕变伸长会超过蠕变伸长限值dLi，结果会处于失效状态．3 算例分析3.1　算例模型上述推导过程阐述了两种极限状态下，索杆结构蠕变伸长限值及安全年限的确定方法．以下通过Matlab软件对一Levy型索穹顶算例进行计算分析，说明具体的求解过程及结果．图1所示索穹顶模型跨度为60 m，设置2道环索，环向12等分．拉索的极限抗拉强度为1 670 MPa，弹性模量为1.95×1011 Pa；压杆的屈服强度为345 MPa，弹性模量为2.06×1011 Pa．恒荷载取0.3 kN/m2，活荷载取0.5 kN/m2．承载能力极限状态下荷载效应应按荷载规范[18]要求取1.3(恒)+1.5(活)进行分析，正常使用极限状态取1.0(恒)+1.0(活)进行分析．各构件的截面积A、钢索直径D、长度L及初始预内力N见表1．10.13245/j.hust.241216.F001图1Levy型索穹顶结构布置示意图10.13245/j.hust.241216.T001表1Levy型索穹顶各类构件基本信息构件标号A/mm2D/mmL/mN/kNRC14632810.19460.02RC24632811.39428.14RC38393811.59711.46DC14032610.44313.97DC24032612.21101.96DC39654011.91486.84HC1691345.18322.76HC22 3406310.721 307.26B14 675—5.00-1 082.62B22 560—5.50-75.18B33 541—6.00-294.49实际工程中的可靠度指标限值通常是基于规范及工程经验进行取值．假定结构的安全等级为一级，按照《建筑结构可靠性设计统一标准》[12]的规定，对于承载能力极限状态，可取β¯=3.7，同时为了对比不同可靠度对于承载能力极限状态下蠕变伸长限值及安全年限的影响，取结构的失效概率为1×10-6，即β¯=4.37进行对比计算．规范中对于正常使用极限状态可取β¯=1.5，同理取β¯=3.5进行对比计算．两种安全系数αc与αd均取1.0．3.2　蠕变伸长限值计算与分析钢索松弛失效模式下各索松弛内力临界值可由式(5)计算得到，即索标RC1～RC3，DC1～DC3，HC1，HC2的松弛内力临界值分别为7.25，8.23，14.90，6.31，7.17，17.10，5.60，39.20 kN；节点变形的临界值取跨度的1/250，即Hck=24 mm进行控制．通过式(21)和(23)可以计算得到变化参数矩阵Mns与Nnn．式(24)所示钢索蠕变模型参数可由钢索松弛失效模式下各索松弛内力临界值及结构整体力学分析得到，两种极限状态部分年数下的蠕变伸长值计算结果见表2和表3，其中最大的年数及蠕变伸长值即对应极限状态的结构安全年限及蠕变伸长限值(以下称为蠕变模型计算结果)．10.13245/j.hust.241216.T002表2各钢索部分年数承载能力极限状态下的蠕变伸长值年限/aRC1RC2RC3DC1DC2DC3HC1HC24815.5917.3719.1415.4817.2220.098.419.944915.7617.5519.3315.6517.4220.288.520.115015.9217.7319.5115.8117.6120.478.620.285116.0817.9119.6915.9817.8020.668.720.465216.2418.0919.8816.1418.0020.848.820.635316.4018.2720.0616.3118.1921.038.920.805416.5618.4520.2416.4718.3821.228.920.975516.7218.6320.4316.6418.5821.419.021.145616.8818.8220.6116.8018.7721.609.121.315717.0419.0020.7916.9718.9621.799.221.485817.2119.1820.9817.1319.1621.989.321.655917.3719.3621.1617.3019.3522.179.421.826017.5319.5421.3417.4619.5422.369.421.996117.6919.7221.5317.6319.7422.559.522.166217.8519.9021.7117.7919.9322.749.622.33mm10.13245/j.hust.241216.T003表3各钢索部分年数正常使用极限状态下的蠕变伸长值年限/aRC1RC2RC3DC1DC2DC3HC1HC23313.7515.2616.8313.4714.3217.277.217.413413.9115.4417.0213.6314.5217.467.317.583514.0715.6217.2013.8014.7117.657.417.753614.2315.8017.3913.9714.9017.847.417.923714.4015.9817.5714.1315.1018.037.518.093814.5616.1617.7514.3015.2918.227.618.263914.7216.3417.9414.4615.4818.417.718.444014.8816.5218.1214.6315.6818.607.818.614115.0416.7018.3114.7915.8718.797.918.784215.2016.8918.4914.9616.0618.988.018.954315.3617.0718.6715.1216.2619.178.019.124415.5217.2518.8615.2916.4519.368.119.294515.6917.4319.0415.4516.6419.558.219.46mm通过Matlab中的fimincon函数求解，得到不同可靠度指标下对应的蠕变伸长限值，对比结果如图2和图3所示．结合图2和表2可以看出，承载能力极限状态下蠕变模型计算结果的安全年限为62 a．相比较而言，β¯=4.37计算出的蠕变伸长限值较为严格，安全年限为49 a；而规范值β¯=3.7时的计算结果则较为宽泛，超过蠕变模型计算结果的蠕变伸长限值及安全年限．若按规范要求维护结构及设定相应的预警机制，则索杆结构会偏不安全．结构在服役期间钢索内力将逐年下降并在特定年数后大概率小于最小预张力，从而突破承载力极限状态．通过计算，当β¯=3.9时的计算结果与数值解基本重合，因此该结构承载能力极限状态的可靠度指标取值应大于3.9．同理，结合图3和表3可以看出，正常使用极限状态下数值解的安全年限为45 a．对比发现β¯=3.5的计算结果相对严格，对应年限为34 a；而规范值β¯=1.5的计算结果较为宽松，超过数值解的蠕变伸长限值及安全年限．若按规范限值控制，则结构在服役期间将大概率发生过大的节点变形，突破结构正常使用极限状态．计算发现，当β¯=3.1时，蠕变伸长限值和安全年限与数值解基本重合，因此该结构正常使用极限状态下的可靠度指标取值应大于3.1．10.13245/j.hust.241216.F002图2承载能力极限状态下不同可靠度指标对应的蠕变伸长限值10.13245/j.hust.241216.F003图3正常使用极限状态下不同可靠度指标对应的蠕变伸长限值4 结论a．基于可靠度理论建立了考虑索杆结构蠕变效应的承载能力极限状态与正常使用极限状态下的极限状态方程和不等式，推导了以长度为变量的内力、位移变化公式，并基于该公式给出了求解极限状态不等式变化参数的具体方法．b．结合钢索蠕变发展模型与数学规划，提出了完整的钢索蠕变伸长限值及安全年限的确定方法．通过算例分析了在结构的两种极限状态下，包含规范取值在内的多种可靠度指标对应的蠕变伸长限值及安全年限，同时与以蠕变发展模型式(24)为基础的计算结果进行了对比，结果表明现有规范取值偏小．因此，在进行此类结构的计算中，可靠度指标应综合考虑规范及对应的结构极限状态取值．