概率密度演化方程差分格式的計(jì)算精度及初值條件改進(jìn)

朱志輝，劉禹兵，高雪萌，周高揚(yáng)，余志武

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

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

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

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

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

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

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

1.1 概率密度演化方程

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

第l組離散代表點(diǎn)集對(duì)應(yīng)的邊界條件和初始條件為：

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

1.2 數(shù)值求解方法

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

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

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

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

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

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

2.1 時(shí)間步長(zhǎng)選取

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

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

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

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

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

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

式中： |Z?|為關(guān)注時(shí)間段的各樣本結(jié)構(gòu)響應(yīng)差值范圍；T為結(jié)構(gòu)響應(yīng)分析總時(shí)長(zhǎng)。

為同時(shí)滿足計(jì)算穩(wěn)定性與精度要求，改進(jìn)時(shí)間步長(zhǎng) dt*應(yīng)為：

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

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

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

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

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

3 差分格式計(jì)算精度的影響因素

3.1 樣本精度影響分析

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

為對(duì)比不同差分格式對(duì)結(jié)構(gòu)可靠度計(jì)算精度的影響，構(gòu)造目標(biāo)三角函數(shù)y(t)：

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

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

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

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

3.2 改進(jìn)初值條件

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

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

圖9(a)中， dt=0.01 s的時(shí)間步長(zhǎng)能夠滿足構(gòu)造樣本精度要求，而對(duì)于 dZ≥1.0075的CFL 條件，過大的空間步長(zhǎng)將導(dǎo)致無法進(jìn)行差分計(jì)算。這進(jìn)一步表明當(dāng)結(jié)構(gòu)響應(yīng)一階時(shí)間導(dǎo)數(shù)過大時(shí)，僅滿足結(jié)構(gòu)響應(yīng)精度要求的時(shí)間步長(zhǎng)并不能滿足概率密度演化方程的差分求解精度。圖9(b)中，初值條件在分布范圍方向上的平滑性很差，這也導(dǎo)致圖9(b)的標(biāo)準(zhǔn)差結(jié)果在初始時(shí)刻大幅震蕩，進(jìn)而導(dǎo)致之后0 s～2 s 內(nèi)的標(biāo)準(zhǔn)差都嚴(yán)重偏離精確解；數(shù)值差分標(biāo)準(zhǔn)差曲線后半段的離散性逐漸增大是因?yàn)橥肯拗破鞯拇嬖趯?dǎo)致概率密度分布范圍只能增大不能減小引起的。

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

以差分解計(jì)算誤差的時(shí)間平均值Erro作為計(jì)算精度指標(biāo)：

式中，Ei為第i時(shí)間點(diǎn)均值或標(biāo)準(zhǔn)差的誤差。

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

4 結(jié)論

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

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

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

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

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