硬化非高斯结构响应首次穿越的Monte Carlo模拟方法

张龙文 卢朝辉 何军 赵衍刚

摘 要：在对比分析已有硬化非高斯模型（Winterstein硬化模型、Ding和Chen模型）的基础上，提出了一个基于Zhao和Lu模型的、新的硬化非高斯模型. 新模型预测偏度和峰度的误差比既有硬化模型小，且最大误差分别为0.311和0.479，表明新模型具有良好精度；同时新模型扩展了Zhao和Lu模型的适用范围. 最后运用新的硬化模型模拟硬化非高斯过程样本，发展了硬化非高斯结构响应首次穿越的Monte Carlo模拟方法. 数值算例验证了本文方法用于硬化非高斯结构响应首次穿越失效概率计算有较高的精度；杭州新火车东站大跨屋盖以及南水北调工程渡槽结构工程实例说明了本文方法的使用过程.

关键词：Winterstein硬化模型；Ding和Chen模型；Zhao和Lu模型；硬化非高斯过程；首穿失效概率

中图分类号：TU311 文献标志码：A

Monte Carlo Simulation Method for the First Passage Probability of Hardening Structural Responses

ZHANG Longwen1， LU Zhaohui1?， HE Jun2， ZHAO Yangang1

（1.College of Civil Engineering， Central South University， Changsha 410075；

2. School of Naval Architecture， Ocean and Civil Engineering， Shanghai Jiaotong University， Shanghai 200240）

Abstract： A new hardening non-Gaussian model based on Zhao and Lu model was proposed， by comparative analysis of the （Wintersteins hardening model， Ding and Chen model）. Skewness error and kurtosis error of the new hardening non-Gaussian model are smaller than the existing hardening response models. The maximum Skewness error and kurtosis error is 0.311， 0.479 respectively. It is indicated that the new hardening non-Gaussian model has good accuracy. At the same time， the new model extends the application range of Zhao and Lu model. Finally， the new hardening non-Gaussian model was applied to simulate hardeing non-Gaussian processes， and Monte Carlo simulation method for the first passage probability of hardening structural responses was developed. Numercial example show that the proposed method has good precision for extimating the first passage probability of hardening structural responses. The long-span roof of Hanzhou New Train Station and the aqueduct of South to North Water Transfer Project were given for illustrating the use process of the proposed method.

Key words： Wintersteins hardening model； Ding and Chen model； Zhao and Lu model； hardening non-Gaussian processes； first passage probability

在结构动力可靠度分析中，首次穿越失效一直是重点研究问题之一. 当工程结构受到诸如地震、风及海浪等隨机荷载激励时，研究结构的首次穿越失效更是具有现实的工程意义与理论价值. 目前，首次穿越失效概率的主要计算方法包括：基于超越率的解析方法[1-2]、基于扩散过程分析的半解析半数值方法[3]、数值积分方法[4]和Monte Carlo模拟方法[5]. 其中Monte Carlo模拟方法为最精确的方法，可以作为“人工试验”来验证其他方法的精度[6-7]. 因此，研究Monte Carlo模拟方法具有十分重要的意义.

Monte Carlo模拟方法的关键是模拟结构的反应样本. 为了节省计算时间，可以从反应的概率分布基础上直接模拟结构反应. 但是，当结构为非线性或荷载为非高斯过程时，只能在结构反应的不完全统计信息基础上进行结构反应过程的模拟. 例如，何军[6]利用结构反应的前四阶矩，基于Winterstein多项式，模拟了非高斯荷载作用下结构的反应过程，并建立了结构首次失效时间分析的模拟方法.

Grigoriu提出的转换过程理论[8]利用转换思想首次将非高斯随机过程变换为标准高斯过程的形式. 此后，发展了多种转换模型用以表示非高斯过程，并在软化反应（具有比高斯分布宽的尾部，即峰度系数>3）取得一定的成果[9]. 对于硬化反应（具有比高斯分布窄的尾部，即峰度系数<3）在实际工程中也经常出现[10-12]，例如，由于海洋波浪或地震作用引起的结构响应、高层建筑围护结构的风压以及风力发电机组的动态响应等，但研究较少. 在已有的转换模型研究中，Winterstein多项式[9]在软化非高斯过程情况运用广泛，但在硬化非高斯过程的运用较少，且许多文献[10-11]发现它的精度不足. 针对Winterstein硬化模型转换精度的不足，基于随机过程的正交展开，Ding和Chen[13]提出了一个更为合理的硬化模型. Ding和Chen对模型的精度进行了量化，但该模型的系数过于复杂. 另外，上述两种模型无法得到类似Winterstein[14]软化非高斯过程模型的形式，不能建立多项式系数与前四阶统计矩（均值、标准差、偏度、峰度）关系的完整表达式. 且当用标准高斯过程表示非高斯过程时，涉及一元三次方程的求根问题. Zhao和Lu[15]提出的四阶矩标准化的函数表达形式，能够建立三次多项式系数与统计矩之间的关系. 该模型简单，避免了上述求根问题，但在硬化非高斯过程的运用中还有待进一步考查. 因此，已有硬化非高斯过程的转换模型还存在不足，其模型的精度及适用范围需进一步调查，以便运用于硬化非高斯结构响应的首穿失效概率计算.

本文基于Zhao和Lu模型的四阶矩标准化函数，提出了新硬化非高斯模型并运用于硬化非高斯过程的首穿失效概率计算. 首先，对已有的Winterstein硬化模型、Ding和Chen模型以及Zhao和Lu模型的精度进行误差分析；接着，基于Zhao和Lu模型，提出了新模型，并进一步讨论了新模型的适用范围；最后，以非线性Duffing振子数值算例验证了本文方法在硬化非高斯结构响应首穿失效概率计算中的有效性；以杭州新火车东站大跨屋盖非高斯脉动风压、南水北调工程渡槽结构的地震反应为例，说明了本文方法的使用过程.

1 硬化非高斯过程的转换模型

1.1 Winterstein硬化模型

Winterstein硬化模型[9]可表达为：

. （1）

式中：U（t）为标准高斯过程；XS（t）=[X（t）-X]/X， h3=3X/6， h4=（4X-3）/24；X， X，3X， 4X 分别为平稳非高斯过程X（t）的均值、标准差、偏度和峰度.

为了对该模型进行误差分析，通过反算硬化非高斯过程的偏度与峰度，并分别与其目标值进行对比分析. 对于非高斯过程X（t）的前四阶中心矩可通过式（2）计算.

记 分别为偏度与峰度的误差，偏度与峰度计算公式表达如下：

， （4a）

. （4b）

根据式（2a）（2b）（3a）（3c）及（4a）（4b）计算模型误差. 图1、图2分别给出了Winterstein硬化模型在3X=0.0、0.2、0.4、0.6不同目标偏度系数下，随着峰度变化的偏度误差曲线和峰度误差变化曲线.

1）图1说明了随着目标偏度的增大，误差变大，且误差曲线的离散性较大.

2）图2说明了在不同目标偏度下，峰度的误差曲线变化趋势基本一致，且随着峰度的增大而减小. 当峰度4X=1.5时，峰度误差达到1.

根据式（1）可知，当3X=-0.2， -0.4， -0.6时，偏度误差及峰度误差的绝对值与3X取正值时相同.

. （6c）

它的应用范围为：

. （7）

根据式（2a）（2b）（4a）（4b）及（6a）（6c）计算模型误差. 图3、图4分别给出了Ding和Chen模型在3X=0.0、0.2、0.4、0.6的不同目标偏度系数下，随着峰度变化的偏度误差曲线和峰度误差曲线.

1）图3说明了该误差曲线离散性较大，且偏度误差的绝对值变化在0～0.75之间.

2）图4说明了在不同目标偏度下，峰度的误差曲线基本一致，且峰度误差的绝对值在0～1之间. 根据等式（5a）（5e）可知，当3X=-0.2、-0.4、-0.6时，偏度误差及峰度误差的绝对值与3X取正值时相同.

1.3 Zhao和Lu模型

Zhao和Lu[15]以一个三次多项式的形式表达：

. （8a）

式中多项式系数l1， k1， k2表达为：

. （8b）

式中l2为：

. （8c）

式（8a）的反函数可表达为：

. （9a）

式中：

； （9b）

. （9c）

根据式（2a）（2b）（4a）（4b）及（8a）（8c）计算模型误差. 图5与图6分别给出了Zhao和Lu模型在3X=0.0、0.2、0.4、0.6的不同目标偏度系数下，随着峰度变化的偏度误差曲线和峰度误差曲线. 图5与图6均表明Zhao和Lu模型在不同目标偏度系数下的误差趋近于0. 因此，通过对比Winterstein硬化模型与Ding和Chen模型的誤差曲线，Zhao和Lu模型显示的误差最小. 根据等式（8a）（8c）可知，当3X=-0.2、-0.4、-0.6时，偏度误差及峰度误差的绝对值与3X取正值时相同. 然而，Zhao和Lu模型的应用范围需满足以下等式：

. （10）

因此，对于在4X≤2.3时的强非高斯硬化过程，该模型将不再适用.

（11b）

式中l2修订为 ，表达为：

. （11c）

图7与图8给出了新模型在3X=0.0、0.2、0.4、0.6的不同目标偏度系数下，随着峰度变化的误差曲线. 图7说明偏度误差绝对值的最大值在0.3左右，图8说明峰度误差绝对值小于0.5. 根据等式（8a）（11a） （11c）及l1可知，当3X=-0.2、-0.4、-0.6时，偏度误差及峰度误差的绝对值与3X取正值时相同.

表1与表2分别给出了Winterstein硬化模型、Ding和Chen模型以及新模型在目标偏度3X=0.0、0.2、0.4、0.6下的偏度误差绝对值最大值、峰度误差绝对值的最大值. 表1说明了本文修正模型偏度误差绝对值的最大值最小，最大误差为0.311. 另外，在3X=0.0时误差为0. Winterstein硬化模型与Ding和Chen模型均有较大的误差，最大误差分别为0.819和0.797.

表2说明了本文修正模型在3X=0.0、0.2、0.4时相比Winterstein硬化模型与Ding和Chen模型的峰度误差绝对值的最大值最小，最大误差为0.479，最小误差为0.428. 而Winterstein硬化模型与Ding和Chen模型的误差绝对值的最大值均有等于1的情况. 因此，根据表1与表2的误差对比分析，说明新模型能提供更高的精度.

2.2 新模型的适用范围

根据式（11c）可知，偏度系数3X与峰度系数a4X的关系应该满足以下等式：

. （12）

一般地，对于常见分布的范围，偏度系数与峰度系数的关系为[16]：

图9显示了新模型、Ding和Chen模型以及式（13）的适用范围. 说明了式（12）涵盖了大部分的硬化非高斯分布的范围，且比较于Ding和Chen模型的适用范围更大.

3 硬化非高斯过程的首穿失效概率

3.1 U-X变换模拟硬化非高斯过程样本

对于标准高斯过程U（t），一般情况下可根据两种

模型生成[17-18]. 基于硬化非高斯过程前四阶统计矩（均值X、标准差X、偏度3X和峰度4X）以及标准高斯过程U（t）、硬化非高斯过程X（t）可以通过式（8a）表达为含有标准高斯过程U（t）的形式. 当采用新模型生成样本时，k1和k2分别为m1和m2. 若采用Winterstein模型以及Ding和Chen模型生成样本，则可分别根据等式（3a）～（3c）以及（6a）～（6c）进行计算. 图10说明了三次多项式生成硬化非高斯过程X（t）样本的过程.

3.2 首穿失效概率

首穿失效概率定义为在[0，T]的时间范围内，随机过程超越界限x至少一次的概率pf（T），表示为[5]：

. （14）

根据U（t）-X（t）变换，可以通过图10方法生成样本X（t），再根据式（15）的分段函数判定样本是否失效. 式（15）表示为：

（15）

首穿失效概率可通过式（16）进行计算：

. （16）

式中：n为Ij的总和；Nsim为结构硬化非高斯过程X（t）的样本数. 图11说明了首穿失效概率的计算过程.

4 算例

4.1 非线性单自由度Duffing振子

考慮一个单边功率谱为1/，受高斯白噪声激励的非线性单自由度的Duffing振子. 它的运动方程表示为：

. （17）

式中：c为阻尼系数；0为自振频率；为控制非线性的参数. 对于该系统的位移反应概率密度函数f（X）有解析解[19]，表达为：

. （18a）

式中：

（18b）

式（18a）说明结构响应X（t）是非高斯的. 根据结构的参数及其位移概率密度函数求解前四阶矩如表3所示（硬化非高斯过程，4X<3）.

根据U-X变换，对应表3的前四阶矩，模拟得到结构的一次反应样本，如图12所示. 考虑结构的界限水平x=2X，利用表3的前四阶矩，并结合上节说明的模拟方法，计算了10 000个样本函数的首穿失效概率. 图13给出了结构在0～50 s的首穿失效概率. 图13同时给出了根据穿越理论计算的解析解[19]结果，以及Ding和Chen模型、Winterstein模型计算结果. 从图13可看出本文方法计算的首穿失效概率能够与解析解计算结果很好地拟合，而运用Ding和Chen模型、Winterstein模型计算结果与解析结果均有较大差异. 图13说明了本文方法计算首超概率的有效性与准确性，并进一步验证了新模型的准确性.

4.2 杭州新火车东站大跨屋盖脉动风压

大跨屋盖结构往往呈现出较强的非高斯特性. Huang等[20]对杭州新火车东站进行了风洞试验研究，并对其大跨屋盖结构的非高斯风压进行数值模拟. 基于杭州新火车东站实验数据，林巍等[21]也研究了大跨度屋盖结构表面风压的非高斯分布特性. 该风洞实验在90°风向角的风压时程的前四阶矩列于表4.

根据U（t）-X（t）变换，对应表4中B44和B46测点的前四阶矩，分别模拟得到风压时程的典型样本如图14（a）（b）所示.

在90°风向角的风压下，测点B44考虑屋盖结构能够承受的界限水平x=-0.75、-0.80两种情况，计算了10 000个样本函数计算的首穿失效概率，图15给出了结构在0～50 s的首穿失效概率. 测点B46考虑屋盖结构能够承受的界限水平x=-0.85、-0.90两种情况，计算了10 000个样本函数的首穿失效概率，图16给出了结构在0～50 s的首穿失效概率.

4.3 渡槽结构的抗震可靠度

南水北调中线工程某渡槽全长114 m，共3跨，每跨38 m. 该渡槽工程位于地震多发带，地震基本烈度为8度. 文献[22]对该实际工程结构进行了结构地震反应分析，并根据数理统计得到了在5 000条人工地震波下（地震波时长为10 s）的绝对加速度反应数据的统计值.

当自振频率0=20 rad/s时，该结构的绝对加速度反应的前四阶矩分别为：X=1.618 661 m/s2，X= 0.064 916 m/s2，X=-0.108 613，X=2.448 231. 结构反应的前四阶矩说明该结构为硬化非高斯过程. 考虑该渡槽的加速度限值为x=、、三种情况，计算得到10 000个反应样本下结构t=5 s的首次穿越概率分别为：0.219、0.109、0.0169. 它们对应的可靠度分别为：0.775 6、1.231 9、2.122 5. 计算结果表明，随着界限值的增大，渡槽结构的可靠度增大显著.

5 结论

1） 基于Zhao和Lu模型，对模型系数进行了修正. 通过对比已有的硬化模型（Winterstein硬化模型、Ding和Chen模型），说明了新模型的偏度、峰度误差最小：最大误差分别为0.311和0.479.

2） 通过数值算例验证了本文新模型运用于硬化非高斯过程样本的模拟以及首穿失效概率计算的有效性与准确性；通过杭州新火车东站大跨屋盖及渡槽结构的实例分析，说明了本文方法在实际工程中 的应用.

3） 本文新模型可应用于实际工程的首穿失效概率计算及工程结构的动力可靠度评估. 另外，由于结构在动力作用下的破坏指标是建立在首次穿越和塑性累积损伤联合效应的基础上的，因此，考虑累积效应的结构动力问题以及首次穿越和累计效应的作用规律需要进一步深入研究.

参考文献

[1] VANMARCK E H. On the distribution of the first passage time for normal stationary process [J]. International Journal of Applied Mechanics， 1975，42（1）： 215-220.

[2] 张颖，易伟建，谭平，等.大震下中间层隔震体系的随机动力可靠性分析[J].湖南大学学报（自然科学版），2009，36（3）：11-15.

ZHANG Y， YI W J， TAN P， et al. The dynamic reliability of mid-story isolation structure under seldomly occurred earthquake [J]. Journal of Hunan University（Natural Sciences）， 2009，36（3）：11-15.（In Chinese）

[3] 甘春标.高斯白噪声激励下拟可积哈密尔顿系统的可靠度[J].振动与冲击， 2006， 25（2）： 145-151.

GAN C B. Study on reliability of quasi integrable Hamiltonian systems under Gaussian white noises [J]. Journal of Vibration and Shock， 2006， 25（2）： 145-151. （In Chinese）

[4] MADSEN P H， KRENK S. An integral equation method for the first passage problem in random vibration [J]. Journal of Applied Mechanics， 1984，51（3）：674-679.

[5] BAYER V， BUCHER C. Importatnt sampling for first passage problems of nonlinear structures [J]. Probabilistic Engineering Mechanics， 1999，14（1/2）：27-32.

[6] 何軍.非高斯荷载作用下结构首次失效时间分析的Monte Carlo模拟方法 [J]. 振动与冲击， 2007， 26（3）：59-60.

HE J. Monte Carlo simulation for first failure time of structures excited by non-Gaussian load [J]. Journal of Vibration and Shock， 2007， 26（3）：59-60. （In Chinese）

[7] 马涌泉，邱洪兴.“顶吸基隔”结构非平稳随机地震反应分析新方法[J]. 湖南大学学报（自然科学版），2015，42（1）：31-39.

MA Y Q， QIU H X. A new method for analyzing the non-stationary random seimic responses of structure with top-absorption and base-isolation [J]. Journal of Hunan University（Natural Sciences）， 2015，42（1）：31-39.（In Chinese）

[8] GRIGORIU M. Crossings of non-Gaussian translation processes [J]. Journal of Engineering Mechanics， 1984， 110（4）： 610-620.

[9] WINTERSTEIN S R. Nonlinear vibration models for extremes and fatigue [J]. Journal of Engineering Mechanics， 1988， 114（10）： 1772-1790.

[10] CHOI M， SWEETMAN B. The Hermite moment model for highly skewen response with application to tension leg platforms [J]. Journal of Offshore Mechanics and Arcitc Engineering， 2010， 132（2）： 950-956.

[11] HUANG M F， LOU W J， CHAN C M， et al. Peak distrubutions and peak factors of wind-induced pressure processes on tall buildings [J]. Journal of Engineering Mechanics， 2013， 139（12）： 1744-1756.

[12] DING J， GONG K， CHEN X. Comparison of statistical extrapolation methods for the evaluation of long-term extreme response of wind turbine [J]. Engineering Structures， 2013， 57（4）： 100-115.

[13] DING J， CHEN X. Moment-based translation model for hardening non-Gaussian response processes [J]. Journal of Engineering Mechanics， 2016， 142（2）：06015006-1-06015006-7.

[14] WINTERSTEIN S R， UDE T C， KLEIVEN G. Springing and slow-drift responses： predicted extremes and fatigue vs. simulation [C]// BOSS-94. Cambrige： Massachusetts Institute of Technology， 1994： 1-15.

[15] ZHAO Y G， LU Z H. Fourth-moment standardization for structural reliability assessment [J]. Journal of Structural Engineering， 2007， 133（7）：916-924.

[16] STUART A， ORD K. Kendalls advanced theory of statistics： distribution theory [M]. 6th ed. New Jersey： Wiley， 2010：22-23.

[17] GRIGORIU M. On the spectral representation method in simulation [J]. Probabilistic Engineering Mechanics， 1993， 8（2）：75-90.

[18] SHINOZUKA M， DEODATIS G. Simulation of stochastic processes by spectral representation [J]. Applied Mechanics Reviews， 1991， 44（4）：191-204.

[19] SOONG T T， GRIGORIU M. Random vibration of mechanical and structural systems [M]. Englewood Cliffs： Prentice-Hall， 1993： 220-285.

[20] HUANG M F， PAN X， LOU W， et al. Hermite extreme value estimation of non-Gaussian wind load process on a long-span roof structure [J]. Journal of Structural Engineering， 2014， 140（9）： 04014061.

[21] 林巍， 黃铭枫， 楼文娟.大跨屋盖脉动风压的非高斯峰值因子计算方法 [J]. 建筑结构， 2013， 43（15）：83-87.

LIN W， HUANG M F， LOU W J. Peak factor method for non-Gaussian fluctuating pressures on long-span roofs [J]. Building Structure， 2013， 43（15）：83-87. （In Chinese）

[22] 易明.地震作用下渡槽结构的动力可靠度分析 [D]. 郑州：郑州大学水利与环境学院， 2010：52-53.

YI M. Seismic dynamic reliablity analysis of the aqueduct structures [D]. Zhengzhou： School of Water Conservancy and Environment， Zhengzhou University， 2010：52-53. （In Chinese）