朱志輝，劉禹兵，高雪萌，周高揚，余志武

(1. 中南大學土木工程學院，湖南，長沙 410075；2. 中南大學高速鐵路建造技術國家工程研究中心，湖南，長沙 410075)

工程結構服役期內(nèi)會受到各類隨機動力作用[1]，結構響應也表現(xiàn)出明顯的隨機性。開展隨機動力作用下的結構可靠度分析對結構性能和安全性評估至關重要。近年來，由李杰和陳建兵[2-3]提出的廣義概率密度演化理論(Generalized probability density evolution theory, GPDET)為結構的隨機響應預測及其可靠度分析提供了新的途徑。

由于實際工程的復雜性，GPDET 的偏微分方程很難獲得精確理論解析解，因此常采用數(shù)值方法求解[4]。這些數(shù)值算法包括有限元法[5]、路徑積分法[6]、等效線性化法[7]、有限差分法[8]等。當采用適當?shù)牟罘指袷綍r，有限差分法在概率密度演化分析中表現(xiàn)出良好性能[9]。目前常用的概率密度演化方程(PDEE)有限差分法按照差分近似微分的處理方式分為單邊差分、L-W 雙邊差分和總變差減小(TVD)性質(zhì)差分[10]三種數(shù)值差分方法。單邊差分和L-W 雙邊差分法分別利用單向差分和中心差分來近似微分結果，兩者分別具有一階和二階計算精度。而TVD 差分法由于引入了通量限制器，避免了L-W 雙邊差分法可能存在的非物理負數(shù)概率解。因此，李杰和陳建兵[10]構造了一類適用于PDEE 求解的TVD 差分格式并將其沿用于其后開展的各項概率密度演化分析。

GPDET 分為概率分配和PDEE 求解兩大關鍵步驟。概率分配對計算精度影響的研究[11-12]已較為成熟。在概率分配滿足計算精度的前提下，YU等[13]基于GPDET 開展鐵路橋梁在隨機荷載條件下的結構響應分析時，通過對比數(shù)論選點法與蒙特卡洛法的激勵樣本計算結果，驗證了車軌橋隨機振動模型的正確性；XU 等[14]基于GPDET 將車輛—軌道相互作用模型與軌道不平順概率模型有效耦合，同樣利用數(shù)論選點法與蒙特卡洛法計算結果對比，驗證了耦合模型的正確性。以往研究中，GPDET 與蒙特卡洛法的對比結果僅體現(xiàn)了概率分配步驟中代表性樣本對樣本空間的表征性，而兩種方法對比的均值和標準差均采用TVD 差分法計算。由于該差分方法計算精度有限，且在局部極值點附近精度降低[15]，這也導致概率密度演化分析中[16-17]，在荷載作用的后半程(卸載階段)，結構響應的概率密度演化標準差結果依舊持續(xù)增大，這樣的計算結果顯然不滿足精度要求，甚至不符合實際情況。因此，PDEE 求解步驟中產(chǎn)生的計算誤差不容忽視，開展基于TVD 差分法的PDEE數(shù)值求解精度分析非常重要。

通常情況下，PDEE 的初值條件為在空間方向可導性較差的脈沖型狄拉克(Dirac)函數(shù)，這將影響到概率密度數(shù)值解和可靠度分析的精確性。為提升TVD 差分方法的求解精度，石晟等[18]提出了余弦型初值條件，余弦型初值雖比脈沖型初值平滑，但空間方向的初始概率賦值范圍受限于余弦函數(shù)周期，且未給出時間差分步長選取原則。

本文針對PDEE，首先對單邊差分、L-W 雙邊差分和TVD 差分三種差分方法開展單樣本和多樣本的差分精度分析；進而提出了TVD 差分的時間步長合理選取方法和正態(tài)分布型初值條件，通過與蒙特卡洛法的數(shù)理統(tǒng)計均值和方差作對比，驗證了本文方法的合理性。

1 概率密度演化方程及數(shù)值求解方法

1.1 概率密度演化方程

結構響應量Z和隨機變量 Θ組成的增廣系統(tǒng)(Z,Θ)是一個保守的隨機系統(tǒng)，隨機性完全來自于Θ，整個過程中，沒有隨機因素的增減。將(Z,Θ)的聯(lián)合概率密度函數(shù)記為pZ,θ(Z,θ,t)，根據(jù)概率守恒原理[19]，對于給定的離散參數(shù)代表點集θl(l=1,2,···,nsel)，PDEE 可寫為一維偏微分方程形式：

第l組離散代表點集對應的邊界條件和初始條件為：

式中： dz為 空間離散步長；nz為空間離散點數(shù)量。

1.2 數(shù)值求解方法

式(1)的偏微分方程很難給出理論解析解，故通常采用數(shù)值差分方法進行求解計算。通過對概率密度值p進 行時間t方向的一階和二階泰勒展開[9]，利用差分方法近似p的一階和二階偏微分，可得：

1.3 差分格式的數(shù)值耗散

由于式(5)和式(6)利用差分代替微分得到線性算子Aγ，由式(17)可知，差分格式、結構響應方向上各時刻概率密度函數(shù)p的可導性、結構響應一階時間導數(shù)Z˙ 、差分求解次數(shù)k(或時間長短)和概率初值條件p0均對數(shù)值差分結果有直接影響。且任一時刻的概率密度值都可由初值條件p0j迭代獲得。

整理得到單邊差分與雙邊差分概率密度演化算子Bo和Bd的表達式分別為：

當 |λZ˙|＜1時，雙邊差分的概率密度演化算子Bd比單邊差分的Bo更接近1。由式(18)可知，只考慮一階泰勒展開的單邊差分法，其數(shù)值耗散大于考慮二階泰勒展開項的雙邊差分法。結合式(14)可知，TVD 差分法的數(shù)值耗散情況介于單邊差分與雙邊差分之間。

2 數(shù)值差分時間步長及正態(tài)分布型初值條件

2.1 時間步長選取

數(shù)值差分的穩(wěn)定性條件為CFL 條件(Courant-Friedrichs-Lewy condition)：

由式(9)和式(19)可得空間離散步長的要求為：

式(20)表明：只要響應離散步長 dZ足夠大，便可保證概率密度演化計算收斂。而當工程結構受地震[20]、軌道不平順[21]和風載荷[22]等劇烈變化的隨機荷載時，則需要更小的時間步長與空間步長來提升概率密度演化結果計算精度。為呈現(xiàn)出概率密度演化結果的峰值變化，并保證結果誤差符合工程要求，必須使各樣本的概率演化曲線在“時間—空間”網(wǎng)格上有明顯的重疊和分離部分，這就對網(wǎng)格離散精度提出了要求。

有限差分法求解過程示意圖如圖1 所示，豎軸代表響應量，樣本A 和樣本B 分別代表樣本空間中某兩條以通量速度為Z˙A和Z˙B在“時間t—空間Z”網(wǎng)格上進行概率密度演化的曲線。

圖1 中，區(qū)域1 的兩樣本差異小于 dZ，兩樣本的概率密度值將在同一空間離散步左右相鄰網(wǎng)格節(jié)點上表征，致使兩樣本結果差異性較低，隨著類似情況的增加，誤差不斷累積，最終的概率密度分布范圍也將比實際情況更寬(標準差結果失真)。同時，由于概率守恒(各時刻概率之和恒為1)，概率密度峰值也將比實際情況更小。因此，只有當空間離散步長 dZ較小時，各樣本的變化差異性才能很好的體現(xiàn)出來(如區(qū)域2)。

使用TVD 差分法求解PDEE 時，由于通量限制器φ取值范圍為[0, 1]，故TVD 差分法的單樣本概率密度演化曲線的離散性只能隨差分迭代次數(shù)的增加逐漸增大或維持離散性不變，這也將導致最終的概率密度演化結果(或標準差結果)誤差較大，甚至是錯誤的。

式中： |Z?|為關注時間段的各樣本結構響應差值范圍；T為結構響應分析總時長。

為同時滿足計算穩(wěn)定性與精度要求，改進時間步長 dt*應為：

式中：T為概率密度演化分析總時長； Ψ為保守系數(shù)，各樣本響應結果差異越小，保守系數(shù)取值越大。

2.2 正態(tài)分布型初值條件

在評價可靠度時，通常將獲得關注量的響應與規(guī)定限值比較，令“W=規(guī)定限值 -關注量響應”。通過假設等價極值事件Wsin(t)，將該虛擬三角函數(shù)的一階時程導數(shù)Wcos(t)代入式(1)，即可得到虛擬隨機過程的概率密度演化結果。當t=π/2 時 ，等價極值事件與實際過程W相等[23]。

由式(23)可知，概率密度pZ的求解精度會直接影響后續(xù)的可靠度分析結果，通常各樣本對應的PDEE 初始條件為式(3)所示的Z0處的脈沖型Dirac 函數(shù)，而Dirac 函數(shù)初值條件在空間方向振蕩明顯，這導致初值條件關于空間方向的可導性對數(shù)值計算結果產(chǎn)生影響[18]。

為進一步提高TVD 差分的計算精度，本文基于正態(tài)分布函數(shù)提出了一種如圖2 所示的初值條件：

3 差分格式計算精度的影響因素

3.1 樣本精度影響分析

3.1.1 單樣本數(shù)值差分精度

為對比不同差分格式對結構可靠度計算精度的影響，構造目標三角函數(shù)y(t)：

式中：A為三角函數(shù)幅值，取值為2；t為時間，選取范圍為0 s～3.15 s。

選取差分時間步長0.0001 s，按式(4)計算單樣本的均值和標準差(該統(tǒng)計特征值無實際物理意義，僅用于判斷演化結果的離散程度)，結果見圖3。算要求(0.01 s)的概率結果離散性明顯大于0.0001 s的標準差結果。當差分步長小于0.004 s 時，TVD差分在第一個周期內(nèi)的標準差均小于5%。按式(22)計算的時間步長為 dt*≤0.05/(5×2×5)=0.001 s，該步長能很好地保證TVD 差分的單樣本計算精度。圖5(b)中，雙邊差分標準差均符合誤差限值要求，這表明相同誤差限值要求下，雙邊差分可采用更大的差分步長，提升計算效率。

3.1.2 多樣本數(shù)值差分精度

圖6 為不同差分時間步長的多樣本統(tǒng)計特征值。圖6(a)中，差分時間步長對多樣本均值不產(chǎn)生影響。圖6(b)中，TVD 差分的標準差結果受差分時間步長影響顯著，0.0008 s 和0.01 s 時間步長的 π/2時刻處標準差值分別為0.098 56 和0.1147，與精確解0.1 的誤差分別為1.44%和14.7%，隨著概率的演化，各差分時間步長的標準差相對于精確解均出現(xiàn)了增大現(xiàn)象，其中 π時刻精確解為0，而0.01 s 步長計算的標準差為0.071 22，計算結果嚴重偏離精確解。圖6(c)中，由于雙邊差分的單樣本的數(shù)值耗散遠小于TVD 差分，故其多樣本標準差結果也更接近于精確解。雙邊差分法以0.0008 s 和0.01 s 時間步長計算的 π/2時刻處標準差值分別為0.099 69 和0.0987，與精確解0.1 的誤差分別為0.31%和1.3%。因此，當計算結果僅關注均值和標準差等統(tǒng)計指標時，使用數(shù)值耗散較小的雙邊差分法，可在較大步長條件下獲得高精度的統(tǒng)計指標，并大幅縮減計算時間。

3.2 改進初值條件

型初值條件零時刻的概率分布較為尖銳，概率密度函數(shù)(Probability density function, PDF)峰值為12.85，而圖7 的精確解初始時刻PDF 峰值為7.97。

選取原始的時間步長為 dt=0.01 s，根據(jù)式(27)和式(28)即可得到100 條目標樣本曲線。這100 條隨機樣本的均值和標準差分別按傳統(tǒng)數(shù)理統(tǒng)計方法和數(shù)值差分法計算，結果見圖9。

圖9(a)中， dt=0.01 s的時間步長能夠滿足構造樣本精度要求，而對于 dZ≥1.0075的CFL 條件，過大的空間步長將導致無法進行差分計算。這進一步表明當結構響應一階時間導數(shù)過大時，僅滿足結構響應精度要求的時間步長并不能滿足概率密度演化方程的差分求解精度。圖9(b)中，初值條件在分布范圍方向上的平滑性很差，這也導致圖9(b)的標準差結果在初始時刻大幅震蕩，進而導致之后0 s～2 s 內(nèi)的標準差都嚴重偏離精確解；數(shù)值差分標準差曲線后半段的離散性逐漸增大是因為通量限制器的存在導致概率密度分布范圍只能增大不能減小引起的。

圖10 正態(tài)分布型初值條件標準差 σ=dZ時的概率密度分布結果。由于采用了正態(tài)分布型初值條件，其初始時刻PDF 峰值為11.78，與圖8 的Dirac 函數(shù)型初值條件相比，PDF 峰值降低，分布范圍方向上的可導性提高，使概率密度分布結果更接近圖7 的精確解。

以差分解計算誤差的時間平均值Erro作為計算精度指標：

式中，Ei為第i時間點均值或標準差的誤差。

使用不同 σ值( σ 以0.2 dZ為間隔)的正態(tài)分布型初值條件開展數(shù)值計算，結果見圖11。圖11 結果表明：改進的正態(tài)分布型初值條件可通過調(diào)節(jié)標準差 σ達到減小均值和標準差計算誤差的目的。隨著 σ的增大，均值累計誤差一直減小；而標準差累計誤差先減小后增大， σ=dZ時達到最小值。與式(27)精確解的初始概率相比，按式(3)Dirac 函數(shù)型初值條件累加求和得到的初始概率分布更集中(即峰值更大，分布范圍更窄)，因此，正態(tài)分布型初值條件在 σ=dZ之前， σ的增大導致初始概率分布的峰值逐漸降低，分布范圍逐漸變寬，初始概率分布向精確解逼近；隨著 σ的進一步增大，雖然初始概率分布的峰值依舊向精確解靠近，但分布范圍卻比精確解更大，導致標準差累計誤差不再減小。

4 結論

本文分析了單邊差分、L-W 雙邊差分和TVD差分三種差分方法在不同差分時間步長條件下的樣本離散性，指出了差分精度的影響因素，提出了TVD 差分法的時間步長合理選取方法和正態(tài)分布型初值條件，通過數(shù)值算例驗證了時間步長選取方法和正態(tài)分布型初值條件的準確性及普適性，并得到如下結論：

(1) L-W 雙邊差分法的單樣本概率演化離散性最小，誤差允許范圍內(nèi)可采用的時間步長最大，當計算結果僅關注均值和標準差等統(tǒng)計指標，而不考慮概率密度分布情況時，建議使用L-W 雙邊差分法開展高效的數(shù)值求解。

(2) 通過對比三種差分方法的數(shù)值耗散情況，得到數(shù)值差分計算精度的影響因素包括：結構響應方向上各時刻概率密度函數(shù)的可導性、結構響應一階時間導數(shù)、差分求解次數(shù) (或時間長度)和概率初值條件。

(3) 基于TVD 差分開展PDEE 數(shù)值求解時，空間離散步長不僅要滿足CFL 條件要求，還應滿足各樣本差異性的最大值限制，本文提出的差分時間步長選取方法有效降低了TVD 差分的單樣本離散性和多樣本概率密度求和導致的累加誤差，并且差分時間步長可根據(jù)樣本差異程度靈活調(diào)整。

(4) 在傳統(tǒng)的脈沖函數(shù)型初值條件基礎上提出了一種更平滑且可以通過標準差 σ控制的正態(tài)分布型初值條件形式。當 σ取值較小時，該分布條件近似于脈沖型Dirac 函數(shù)初值分布。合適的標準差σ可以提高正態(tài)分布型初值條件的數(shù)值差分精度。