基于可靠度的木结构荷载持续作用效应

何敏娟，张婷钰，倪 春，李 征，陶 铎

（1.同济大学土木工程学院，上海200092；2.加拿大林产品创新研究院，温哥华V6T1Z4）

木材在长期荷载作用下，随着时间的增加材料强度逐渐降低，最终发生破坏。这种随着时间的增加强度下降的现象叫做荷载持续作用效应（duration of load，DOL）。早在1741年，法国海军工程师Buffon就注意到了这个现象，他认为木梁的抗弯承载力不应超过短期强度的50%。20世纪40年代，Wood用花旗松进行了清材小试件的受弯试验，提出了著名的Madison曲线。直到今天，Madison曲线也还在被用来预测木材在DOL效应下的强度退化情况。由于通过长期试验测出木材的长期强度较为费时费力，需要通过数学模型来对长期强度作出合理的预测。为了预测木材的DOL效应，研究者们提出了大量模型，包括经验模型（empirical model）、累积损伤模型（damage accumulation model）、断裂力学模型［1］（fracture mechanics model）、变形动力学模型［2］（deformation dynamic model）和能量模型［3］（energy model）来表示长期强度与短期强度之比和荷载持续时间的关系。1979—1987年，美国学者Gerhards提出了指数型累积损伤模型，先利用清材小试件对模型进行校准，随后又用了足尺试件的长期试验数据进行校准。1978年，加拿大学者Foschi等在Gerhards模型基础上，引入了二阶损伤变量及应力阈值σ0。他们基于以恒定速率增加到某一数值再保持不变的荷载模式对清材小试件进行试验，并开发了BF损伤模型。随后将模型推广到更复杂的荷载模式，并用足尺试件的长期荷载数据进行了模型校准。1986年Foschi和Yao［4］在BF模型的基础上作了改进，提出Foschi-Yao模型，用应力值表示损伤累积，模型参数增加到4个。2017年，中国学者祝恩淳教授曾对木材可靠度进行分析［5］，通过对木构件的γR进行研究进而确定木材强度设计值。本文采用Foschi-Yao累积损伤模型，基于我国恒荷载、楼面活荷载及雪荷载的概率统计模型，对荷载持续作用效应进行研究。

1 在短期荷载作用下的可靠度计算

1.1 功能函数

可靠度分析就是要根据可靠性总原则的规定和我国可靠度分析的一般要求，结合木材的强度特性，建立极限状态方程及功能函数。

恒荷载和活荷载作用下的荷载效应组合E（Q+D）可以表示为

式中：d为归一化恒荷载，d=D/Dn，D为恒荷载真实值，Dn为恒荷载标准值；q为归一化活荷载，q=Q/Qn，Q为活荷载真实值，Qn为活荷载标准值；ρ=Qn/Dn，ρ一般取值为0.5、1.0、1.5、2.0、3.0、4.0。

中国规范《木结构设计标准：GB50005—2017》［6］中的极限状态设计公式为

式中：γG和γQ分别为永久荷载和可变荷载的分项系数；E（Dn）、E（Qn）分别为永久荷载和可变荷载的荷载效应；f（0.05）为材料强度特征值，即fk，短期强度分布的5%分位值；KDOL为长期荷载影响系数，在不考虑长期荷载效应时取定值1.0。

结合式（1）和式（2）可以得到

则功能函数G可以表示为

构件抗力R在可靠度计算中为随机变量，假定服从对数正态分布。根据GB50005的规定，R为

式中：长期荷载影响系数KQ3=KDOL，在不考虑长期荷载作用时取值1.0；构件几何特征不定性影响系数KA、计算模式误差影响系数KP及木材或木产品的短期强度f均为随机变量。f为足尺试验强度，服从对数正态分布。抗力R可根据随机变量理论，得到其均值和变异系数。

中国规范《木结构设计标准：GB50005—2017》中的极限状态设计公式为

1.2 荷载分布及参数

在进行短期荷载作用下的可靠度计算时，只需考虑荷载变量的分布类型及其相应的归一化均值和变异系数，不考虑荷载的时变效应。其中归一化均值为荷载平均值与标准值之比。根据《木结构设计标准：GB50005—2017》中的规定，各项荷载统计参数及分布类型如表1所示。

表1 荷载统计分布及其参数Tab.1 Distribution type of loads and corresponding parameters

表1中，我国规范对雪荷载的统计分布参数的规定忽略了城市间的差异，均采用统一的分布参数值。但为了保证短期雪荷载模型与长期雪荷载模型间的统一，在进行短期荷载作用下的可靠度计算时，不同城市的雪荷载（50年重现期）的统计参数均采用相应的长期雪荷载模型，计算得到的结果如表2。

表2 典型城市的雪荷载分布参数Tab.2 Parameters of snow load distribution in typi⁃cal cities

1.3 可靠度分析

研究采用国际安全度联合委员会（The Joint Committee of Structural Safety，JCSS）推荐使用的改进的一次二阶矩法，即JC法，计算可靠度指标β与γR之间的关系［7］。根据文献［8］的规定，针对不同安全等级，目标可靠度指标按表3分别取为β0=2.7、3.2、3.7、4.2。本研究在安全等级为二级状态下进行。

表3 结构构件承载力极限状态下的可靠指标Tab.3 Reliability index of the structure limit state

此外，在建立γR与材料强度变异系数COV的关系时，需要利用满足对数正态分布的材料抗力均值与五分位值的关系。

将功能函数重新表示为

该种木产品期望达到的目标可靠度指标以及抗力、恒荷载和活荷载等各随机变量的均值和变异系数均已知，通过迭代可以计算出功能函数中项的均值代入材料抗力均值与五分位的关系式（8），可以得到

进一步可以得到γR与COV的关系

2 在长期荷载作用下的可靠度计算

2.1 损伤模型

累积损伤模型通过引入并没有实际物理意义的损伤变量，在一定荷载历史条件下，对结构单元的累积损伤进行定量计算，从而预测单元的失效时间。常用的共有3种损伤累积模型［9］：Gerhards模型、Barrett-Foschi模型、Foschi-Yao模型，本文采用Foschi-Yao模型，即加拿大模型［10］进行研究，对应力及强度值进行归一化，即统一除以材料强度特征值，并采用以应力比为自变量的函数形式。

式中：α为损伤状态变量；τ（t）为归一化后任意时刻的应力。a、b、c、n和σ0这5个模型参数假定为对数正态分布，对于一个给定的结构构件而言为常数，但对于不同的构件而言为随机变量。τs为构件归一化后的短期标准强度，呈对数正态分布，材料的实际短期标准强度通过加载速率为ks的短期斜坡加载试验测得，同样对于一个给定的结构构件而言为常数，但对于不同的构件而言为随机变量。模型参数σ0定义为界限应力比，实际应力比超过界限应力比时材料才会发生损伤累积，即当-σ0≤0时，材料不会发生损伤累积。

2.2 损伤计算

功能函数G可以表示为

初始状态时损伤变量α=0，即功能函数G=1；当损伤变量α=1，即功能函数G=0时，材料发生破坏。

对于任意荷载时程，一般无法求得式（12）的积分表达式，也无法直接计算破坏时间Tf（即该时刻损伤状态变量α为1）并与使用年限进行对比来判断材料是否会破坏。因此采用近似的方法，将荷载时程细分为时间间隔为Δt的小段（在实际分析过程中，根据活荷载、雪荷载变化的时间间隔分段，每段时长可能不同，但段内荷载保持不变），设荷载时程中第i段的初始损伤为αi-1，则该段荷载结束时的损伤αi可表示为

其中

参数a可写作变量ks、τs、b和σ0的函数，φi为荷载作用下的材料应力比，具体表示为

式中：变量d、q、ρ、γR、γG、γQ和式（4）中的含义相同。荷载时程中任意时刻的损伤均能通过式（14）的递推关系求得。

2.3 荷载模型

2.3.1 恒荷载

恒荷载设计值一般根据材料的平均重量计算。为了考虑结构在使用年限内重量变化的不确定性，一般假定恒荷载D是正态随机变量。恒荷载与其标准值的比值d可表示为

2.3.2 楼面活荷载

在我国荷载规范中，设计基准期内住宅楼面活荷载的均值要明显高于办公室的楼面活荷载均值，同时住宅楼面活荷载的变异系数要更小，因此本文中的荷载模型以住宅的统计参数为依据。

我国《建筑结构荷载规范：GB50009—2012》［11］规定，楼面活荷载按其随时间变异的特点，可分持续性和临时性两部分。出于分析上的方便，对2种活荷载的分布类型采用极值I型分布，参考国外的研究方法，采用JCSS（结构安全联合委员会）提出的楼面活荷载模型拟合实际情况［12］。

持续性活荷载和临时性活荷载模型分别采用以下假定：

（1）持续性荷载出现的时刻X1，X2，…服从泊松分布，所以荷载变化的时间间隔服从均值为λsus的指数分布。

（2）持续性荷载的幅值Psus服从均值为μsus、标准差为σsus=的伽马分布，其中A0为FBC-field（Ferry Borges-Castanheta模型域）的参考面积，A为从属面积，κ为形状系数。

（3）临时性荷载出现的时刻X1，X2，…也服从泊松分布，所以荷载变化的时间间隔服从均值为λtem的指数分布。

（4）临时性荷载的幅值Ptem服从均值为μtem、标准差为σtem=的伽马分布，临时性荷载的时长为ttem。

为了得到基准期为50年的最大组合荷载的分布，将进行5 000次重复的计算机模拟，每次模拟采用不同的持续性和临时性荷载的序列组合。通过2种荷载序列的叠加可以计算每次模拟下的最大组合荷载。选取模拟数据中较大的50%，利用极值I型分布拟合最大荷载值的模拟结果，得到活荷载的相关分布参数。

2.3.3 雪荷载

年最大场地雪荷载G服从极值I型分布，其中，分布参数可以用平均值和变异系数VG得到。则服从极值I型分布的随机变量g为

参数p为0到1间的均匀分布随机变量，其他参数可以表示为

式中：A、B为与位置参数μ和尺度参数α有关的参数，可以根据不同城市10年、50年和100年重现期下的基本雪压值计算得出。由于我国不同城市的基本雪压差距较大，需对各个城市分别进行计算，本文以南京、哈尔滨、北京、上海、金华为例进行说明。假定一年有Ns段等时长的荷载段，每一段的荷载保持不变且相互独立。为了概率分布的统一，该分布必须满足Ns个荷载段中荷载最大值的分布与年荷载最大值的分布相同。通过类比可以得到

式中：pe为一个荷载段下雪的概率。

2.4 可靠度分析

采用蒙特卡洛（Monte Carlo）方法［13］计算考虑长期荷载效应下构件的可靠度。为确保计算结果精确性，样本量取值为NR=100 000次。恒荷载与活荷载组合的计算方法如下：

（1）选择γR、活荷载标准值与恒荷载标准值之比ρ以及材料强度变异系数。

（2）为每次重复模拟建立载荷序列。①选择一个标准正态随机数来计算恒荷载d，并假定在设计周期内保持恒定。②生成活荷载。对于楼面活荷载：反复生成2个均匀分布随机数来计算持续性荷载qsus的大小和持续时间，直到总的持续时间等于或大于使用年限。反复生成3个均匀分布随机数来计算临时性荷载qtem的大小、持续时间及时间间隔，直到总的时间等于或大于使用年限。对于雪荷载：对于每一个荷载段，生成一个0到1之间的均匀分布随机数。如果该随机数大于pe，说明该荷载段无雪荷载；反之，则再生成一个0到1之间的均匀分布随机数，用于计算雪荷载q。

（3）计算每个样本在设计周期内的损伤累积。①根据优化得到的参数b、c、n、σ0和τs的对数正态分布特征选取5个对数正态分布随机数作为相应的取值。②根据递推关系式（14）计算累积损伤变量α。③计算功能函数

如果G≥0，说明样本未破坏；如果G＜0，则说明样本发生破坏。

（4）接下来将上述步骤（2）、（3）重复进行NR次。

（5）统计在使用年限i年内样本破坏的数量NFi并计算破坏概率Pf（i）。

（6）计算相应的可靠度指标βi。

（7）根据要求，可以改变γR、ρ以及COV的取值并从步骤（1）开始进行重复计算。

3 确定损伤模型参数

加拿大林产品创新研究院（FP Innovations）对Hemlock进行了5组长期试验，持荷时间分别为3个月、1年和4年，长期荷载值为20.68 Mpa和31.03 Mpa。试验时试件的应力时程如图1所示。

试验时较为重要的是试件的破坏时间Tf（即该时刻损伤状态变量α为1）。加载时，荷载线性增加到一个定值τc，然后保持不变直到试件破坏（如图1），通过式（12）的积分可得破坏时间的表达式为

式中：tc为试验开始长期持荷时所对应的时间；αc为此时的累积损伤。

图1 试验的应力时程Fig.1 Specific stress history used in tests

由于长期试验的结果重复性较好，取其中试验时间较长的2组数据作为依据对损伤模型进行拟牛顿法非线性优化（nonlinear optimization）［14］。优化的过程大致如下，由MATLAB软件编程实现。

（1）选定各参数均值和变异系数的初值，一般来说初值的数量级不能跟最优解的数量级相差太多，选用的初值如表4。

表4 各参数的初值选择Tab.4 Initial values chosen for model parameters

（2）根据变量b、c、n、σ0和τs的分布生成随机数序列，重复2 000次然后代入式（26）计算破坏时间Tf，按照大小进行排序得到相应的累积概率分布。

（3）将模拟得到的累积分布与试验数据进行对比，其差值由函数Ф表示

式中：N为用于比较时选取的概率水平的数量；Tfi、Tdi分别为相同概率水平下的模拟得到的破坏时间和试验得到的破坏时间。

（4）保持变量随机数的次序不变，分别改变不同变量的均值和变异系数并重复步骤（2）、（3），求出式（27）中函数Ф相对于各变量均值和变异系数的偏导函数并用于接下来的优化程序。

（5）利用自编的非线性优化算法寻找式（27）的最小值，即可得到对应的各参数均值和变异系数的最优解。

表5 各参数的优化结果Tab.5 Optimization results of model parameters

通过优化得到损伤累积模型的参数如表5所示，图2及图3对比了由该损伤模型计算得到的和长期试验得到2种荷载水平下的累积破坏概率-对数破坏时间关系。从图中可以看出，该损伤模型的拟合效果较好，而且该模型能够有效预测持荷时间超过4年后的累积破坏概率增长情况。

图2 荷载为31.03MPa时试验和模拟的累积破坏概率-对数破坏时间对比Fig.2 Damage accumulation model versus experi⁃mental results at 31.03MPa

图3 荷载为20.68 Mpa时试验和模拟的累积破坏概率-对数破坏时间对比Fig.3 Damage accumulation model versus experi⁃mental results at 20.68 Mpa

4 荷载持续时间影响系数

计算的基本过程如图4所示，该图表明2种情况下γR和β之间的关系：短期荷载曲线（曲线1）不考虑材料损伤累积，只考虑材料的短期强度，由第1节方法得到；长期荷载曲线（曲线2）则考虑荷载的持续时间影响，由第2节可靠度分析的方法得到。对于相同的目标可靠度，曲线1（不考虑材料损伤累积）对应的抗力分项系数为γR，1，而曲线2（考虑材料损伤累积）对应的抗力分项系数为γR，2，则可按式（28）计算荷载持续时间影响系数KDOL。

图4 计算荷载持续时间影响系数示意Fig.4 Procedure to determine the DOL adjustment factor K DOL

4.1 恒荷载

对于只有恒荷载的工况，在模拟中将活荷载与恒荷载的比率ρ设为0.001，对COV由0.05～0.50变化进行全面分析。对于受弯、受压以及受拉构件，在β=3.2、β=3.7的情况下，γR，1、γR，2与COV的关系如图5所示。受弯及受压构件的结果相近，在β=3.2的条件下，当COV由0.05增加到0.50时，γR，1由1.026增加到1.689，γR，2由2.231增加到2.856，根据式（28）计算可得KDOL值由0.46增加到0.59左右。对于受拉构件，在β=3.7的条件下，当COV由0.05增加到0.50时，γR，1由1.082增加到2.148，γR，2由2.489增加到3.707，根据式（28）计算可得KDOL值由0.43增加到0.58左右。可以发现，恒荷载工况下的KDOL值随着COV的增大而增大，变化范围在0.43到0.59之间，综合后续工况结果，本文取值为0.52。

4.2 恒荷载且楼面活荷载

图5 ρ=0.001时抗力分项系数与材料强度变异系数的关系Fig.5 γR versus C OV atρ=0.001

对于恒荷载且楼面活荷载的工况，在模拟中分别考虑楼面活荷载与恒荷载的比率ρ为0.5、1.0、1.5、2.0、3.0和4.0时的情况，对0.05～0.50范围内的变化进行全面分析。其中ρ=1.0时的γR，1、γR，2与COV关系如图6，受弯及受压构件的结果相近，在β=3.2的条件下，当COV由0.05增加到0.50时，γR，1由0.958增加到1.412，γR，2由1.493增加到2.036，根据式（28）计算可得KDOL值由0.64增加到0.69左右。对于受拉构件，在β=3.7的条件下，当COV由0.05增加到0.50时，γR，1由1.048增加到1.802，γR，2由1.625增加到2.617，根据式（28）计算可得KDOL值也由0.64增加到0.69左右。对其他ρ进行同样计算分析后，可以认为在恒荷载且楼面活荷载的工况下，对于受弯、受压及受拉构件，在β=3.2和β=3.7情况下的KDOL值基本相同，且KDOL随着ρ的增大而增大，直到ρ大于3.0时，KDOL基本保持在0.80左右。在ρ=0.5、ρ=1.0的情况下，KDOL也随着COV的增大而增大。在ρ=1.5、ρ=2.0、ρ=3.0、ρ=4.0的情况下，KDOL与COV无明显关系，COV在0.05～0.50范围内变化时，KDOL的变化很小。可以得出，COV对KDOL的影响随着ρ的增大而明显减小。当ρ＞1.0时，可以认为这时COV对KDOL几乎没有影响。

图6 ρ=1.0时抗力分项系数与材料强度变异系数的关系Fig.6 γR versus C OV atρ=1.0

鉴于只有恒荷载的工况等同于荷载比率ρ为0.001时的恒荷载且活荷载工况，综合考虑可近似得

4.3 恒荷载且雪荷载（南京）

对于恒荷载且雪荷载的工况，在模拟中考虑城市南京的雪荷载与恒荷载比率ρ为0.5、1.0、1.5、2.0、3.0和4.0等情况，并对材料强度变异系数在0.05～0.50范围内变化进行全面分析。其中ρ=4.0时的γR，1、γR，2与COV关系如图7，受弯、受压及受拉构件的结果相近，在β=3.2和β=3.7的条件下，由式（28）可得，当COV在0.05到0.50范围内变化时，KDOL值在0.79到0.81范围内波动。对其他活、恒荷载比率进行同样计算分析后，可以认为在恒荷载且雪荷载（南京）的工况下，受弯、受压及受拉构件结果相近。当ρ大于1.0时，KDOL基本保持在0.80左右。同时，COV对KDOL的影响随着ρ的增大而明显减小。当ρ大于0.5，COV在0.05～0.50范围内变化时，KDOL的变化很小，可认为这时COV对KDOL几乎无影响。

图7 ρ=4.0时南京抗力分项系数与材料强度变异系数的关系Fig.7 γR versus C OV atρ=4.0 in Nanjing

鉴于只有恒荷载的工况亦等同于荷载比率ρ为0.001时的恒荷载且雪荷载工况，综合考虑可以近似得到南京的KDOL，NJ与ρ的关系式如式（30）所示。

4.4 恒荷载且雪荷载（其他城市）

对于恒荷载且雪荷载的工况，在模拟中考虑城市哈尔滨、北京、上海、金华的ρ为0.5、1.0、2.0和4.0等情况，并对材料强度变异系数由0.05～0.50变化进行全面分析。其中哈尔滨地区ρ=4.0时的γR，1、γR，2与COV关系如图8所示，受弯、受压及受拉构件的结果相近，在β=3.2和β=3.7的条件下，根据式（28）计算可得，当COV在0.05到0.50范围内变化时，KDOL值在0.80到0.82范围内波动。对其他活、恒荷载比率及其他城市进行同样计算分析后，可以认为在恒荷载且雪荷载（哈尔滨、北京、上海、金华）的工况下，受弯、受压及受拉构件结果相近。当ρ大于1.0时，KDOL基本保持在0.81左右。同时，COV对KDOL的影响随着ρ的增大而明显减小。当ρ大于0.5，COV在0.05～0.50范围内变化时，KDOL的变化很小，可以认为这时COV对KDOL几乎没有影响。

图8 ρ=4.0时哈尔滨抗力分项系数与材料强度变异系数的关系Fig.8 γR versus C OV atρ=4.0 in Harbin

因只有恒荷载的工况等同于ρ=0.001时恒荷载且雪荷载工况，综合考虑可近似得哈尔滨（HEB）、北京（BJ）、上海（SH）、金华（JH）的KDOL与ρ关系

5 结论

根据中国规范GB50005—2017及GB50009—2012中荷载的概率统计模型，对荷载的时变过程进行模拟，采用Foschi-Yao模型，将短期及长期荷载下的可靠度计算结果进行对比，得到不同条件下的荷载持续作用效应的影响情况。具体结论如下：

（1）在恒荷载工况下，KDOL值由0.430增加到0.590，综合后续工况结果，本文取值0.520。

（2）在恒荷载且住宅楼面活荷载工况下，ρ＜3.00时KDOL值为0.670+0.280lgρ，ρ＞3.00时KDOL值为0.800。

（3）在恒荷载且雪荷载工况下，南京地区ρ＜1.50时KDOL值为0.745+0.370lgρ，ρ＞1.50时KDOL值为0.810；哈尔滨地区ρ＜2.00时KDOL值为0.745+0.230lgρ，ρ＞2.00时KDOL值为0.810；北京地区ρ＜2.00时KDOL值为0.745+0.250lgρ，ρ＞2.00时KDOL值为0.820；上海地区ρ＜1.00时KDOL值为0.800+0.330lgρ，ρ＞1.00时KDOL值为0.800；金华地区ρ＜2.00时KDOL值为0.760+0.200lgρ，ρ＞2.00时KDOL值为0.820。

作者贡献申明：

何敏娟：提供目前国内外研究进展情况，把控研究整体思路方向，审查研究内容及结果的合理性。

张婷钰：负责具体的计算分析工作，整理论文思路并撰写论文。

倪春：提供研究背景及主要研究思路方法，审查研究内容及结果的合理性。

李征：提供研究思路，审查论文内容，把控研究方向。

陶铎：负责将模型抽象化原理转换为Matlab语言，为后续计算分析做铺垫。