張龍文 盧朝輝 何軍 趙衍剛

摘 要：在對比分析已有硬化非高斯模型（Winterstein硬化模型、Ding和Chen模型）的基礎上，提出了一個基于Zhao和Lu模型的、新的硬化非高斯模型. 新模型預測偏度和峰度的誤差比既有硬化模型小，且最大誤差分別為0.311和0.479，表明新模型具有良好精度；同時新模型擴展了Zhao和Lu模型的適用范圍. 最后運用新的硬化模型模擬硬化非高斯過程樣本，發(fā)展了硬化非高斯結構響應首次穿越的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模擬方法的關鍵是模擬結構的反應樣本. 為了節(jié)省計算時間，可以從反應的概率分布基礎上直接模擬結構反應. 但是，當結構為非線性或荷載為非高斯過程時，只能在結構反應的不完全統(tǒng)計信息基礎上進行結構反應過程的模擬. 例如，何軍[6]利用結構反應的前四階矩，基于Winterstein多項式，模擬了非高斯荷載作用下結構的反應過程，并建立了結構首次失效時間分析的模擬方法.

Grigoriu提出的轉換過程理論[8]利用轉換思想首次將非高斯隨機過程變換為標準高斯過程的形式. 此后，發(fā)展了多種轉換模型用以表示非高斯過程，并在軟化反應（具有比高斯分布寬的尾部，即峰度系數>3）取得一定的成果[9]. 對于硬化反應（具有比高斯分布窄的尾部，即峰度系數<3）在實際工程中也經常出現[10-12]，例如，由于海洋波浪或地震作用引起的結構響應、高層建筑圍護結構的風壓以及風力發(fā)電機組的動態(tài)響應等，但研究較少. 在已有的轉換模型研究中，Winterstein多項式[9]在軟化非高斯過程情況運用廣泛，但在硬化非高斯過程的運用較少，且許多文獻[10-11]發(fā)現它的精度不足. 針對Winterstein硬化模型轉換精度的不足，基于隨機過程的正交展開，Ding和Chen[13]提出了一個更為合理的硬化模型. Ding和Chen對模型的精度進行了量化，但該模型的系數過于復雜. 另外，上述兩種模型無法得到類似Winterstein[14]軟化非高斯過程模型的形式，不能建立多項式系數與前四階統(tǒng)計矩（均值、標準差、偏度、峰度）關系的完整表達式. 且當用標準高斯過程表示非高斯過程時，涉及一元三次方程的求根問題. Zhao和Lu[15]提出的四階矩標準化的函數表達形式，能夠建立三次多項式系數與統(tǒng)計矩之間的關系. 該模型簡單，避免了上述求根問題，但在硬化非高斯過程的運用中還有待進一步考查. 因此，已有硬化非高斯過程的轉換模型還存在不足，其模型的精度及適用范圍需進一步調查，以便運用于硬化非高斯結構響應的首穿失效概率計算.

本文基于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 分別為平穩(wěn)非高斯過程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）

式中l(wèi)2為：

. （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）

式中l(wèi)2修訂為 ，表達為：

. （11c）

圖7與圖8給出了新模型在3X=0.0、0.2、0.4、0.6的不同目標偏度系數下，隨著峰度變化的誤差曲線. 圖7說明偏度誤差絕對值的最大值在0.3左右，圖8說明峰度誤差絕對值小于0.5. 根據等式（8a）（11a） （11c）及l(fā)1可知，當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]. 基于硬化非高斯過程前四階統(tǒng)計矩（均值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為自振頻率；為控制非線性的參數. 對于該系統(tǒng)的位移反應概率密度函數f（X）有解析解[19]，表達為：

. （18a）

式中：

（18b）

式（18a）說明結構響應X（t）是非高斯的. 根據結構的參數及其位移概率密度函數求解前四階矩如表3所示（硬化非高斯過程，4X<3）.

根據U-X變換，對應表3的前四階矩，模擬得到結構的一次反應樣本，如圖12所示. 考慮結構的界限水平x=2X，利用表3的前四階矩，并結合上節(jié)說明的模擬方法，計算了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. 該渡槽工程位于地震多發(fā)帶，地震基本烈度為8度. 文獻[22]對該實際工程結構進行了結構地震反應分析，并根據數理統(tǒng)計得到了在5 000條人工地震波下（地震波時長為10 s）的絕對加速度反應數據的統(tǒng)計值.

當自振頻率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） 本文新模型可應用于實際工程的首穿失效概率計算及工程結構的動力可靠度評估. 另外，由于結構在動力作用下的破壞指標是建立在首次穿越和塑性累積損傷聯(lián)合效應的基礎上的，因此，考慮累積效應的結構動力問題以及首次穿越和累計效應的作用規(guī)律需要進一步深入研究.

參考文獻

[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] 甘春標.高斯白噪聲激勵下擬可積哈密爾頓系統(tǒng)的可靠度[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] 馬涌泉，邱洪興.“頂吸基隔”結構非平穩(wěn)隨機地震反應分析新方法[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]. 鄭州：鄭州大學水利與環(huán)境學院， 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）