劉雪瑩, 譚曉慧, 費(fèi)鎖柱, 侯曉亮

(合肥工業(yè)大學(xué) 資源與環(huán)境工程學(xué)院,安徽 合肥 230009)

在進(jìn)行巖土工程的隨機(jī)有限元分析時(shí),通常利用隨機(jī)場來表示巖土體參數(shù)的空間變異性。一個(gè)隨機(jī)場H(x)就是一個(gè)隨機(jī)過程,其自變量是連續(xù)空間位置參數(shù)x。給定一個(gè)空間位置x0,H(x0)為一個(gè)隨機(jī)變量。隨機(jī)場的實(shí)現(xiàn)就是一個(gè)確定的連續(xù)函數(shù)。隨機(jī)場模型一般不能在隨機(jī)有限元分析中直接使用,應(yīng)先將隨機(jī)場用有限個(gè)隨機(jī)變量表示,這個(gè)過程稱為隨機(jī)場的離散。常用的隨機(jī)場離散方法有中心點(diǎn)法、局部平均法、Karhunen-Loéve級數(shù)展開法(簡稱“KL展開法”)等[1-3]。其中,采用KL展開法進(jìn)行隨機(jī)場的離散能得到連續(xù)的隨機(jī)場(即隨機(jī)場的離散值在網(wǎng)格邊界不會(huì)發(fā)生突變),而且具有離散效率高、得到的隨機(jī)變量數(shù)較少等優(yōu)點(diǎn)。因此,KL展開法是一種較為理想的隨機(jī)場離散方法。

在采用KL展開法進(jìn)行隨機(jī)場離散時(shí),必須確定級數(shù)的展開項(xiàng)數(shù)(又稱“截?cái)囗?xiàng)數(shù)”),KL展開法的展開項(xiàng)數(shù)即離散后的隨機(jī)變量數(shù)。理論上,展開項(xiàng)數(shù)越多,離散結(jié)果越精確,但是,由于隨機(jī)場離散時(shí)計(jì)算量也隨著展開項(xiàng)數(shù)的增加而增加,以及后續(xù)的可靠度分析計(jì)算量也隨之增大,因此,必須合理確定級數(shù)的展開項(xiàng)數(shù)。文獻(xiàn)[4-6]采用隨機(jī)場期望能的比率因子來確定展開項(xiàng)數(shù);文獻(xiàn)[7]采用逐點(diǎn)均方誤差的平均值來確定離散誤差及展開項(xiàng)數(shù)。本文通過理論推導(dǎo)及實(shí)例計(jì)算,分析了離散誤差的影響因素,探究了這2類誤差計(jì)算公式之間的關(guān)系,提出了KL展開法誤差計(jì)算的新公式。

1 隨機(jī)場離散的KL展開法

1.1 KL展開法的基本原理

記平穩(wěn)正態(tài)隨機(jī)場的均值與均方差分別為μ、σ,自相關(guān)函數(shù)為ρ(x1,x2),則隨機(jī)場H(x)可展開為:

(1)

其中,ξi為M個(gè)互為獨(dú)立的標(biāo)準(zhǔn)正態(tài)變量;M為隨機(jī)場的展開項(xiàng)數(shù)(截?cái)囗?xiàng)數(shù));λi為自相關(guān)函數(shù)ρ(x1,x2)前M個(gè)由大到小排列的特征值;fi(x)為與λi對應(yīng)的特征函數(shù)。

將fi(x)用正交基函數(shù)hk(x)展開[1,7],即

(2)

其中,dik為特征向量(即特征向量矩陣D的第i列向量)的第k個(gè)元素。經(jīng)變換,λi、dik計(jì)算公式為:

AD=DΛ

(3)

其中,Λ為M個(gè)特征值組成的對角線矩陣;A為M×M階方陣,其元素計(jì)算公式為:

(4)

其中,j為第n個(gè)離散點(diǎn);hk(x)為fi(x)的正交基函數(shù)展開項(xiàng),可用Legendre多項(xiàng)式函數(shù)表示。

k=1,2,…,∞

(5)

其中,Lk-1(·)為標(biāo)準(zhǔn)Legendre多項(xiàng)式函數(shù);T、a分別為平移參數(shù)與縮放參數(shù),T=(xmax+xmin)/2,a=(xmax-xmin)/2。此時(shí),(4)式可簡寫為:

(6)

對于空間二維問題,特征值、特征函數(shù)可分別由x、y2個(gè)坐標(biāo)軸方向的特征值乘積、特征函數(shù)乘積表示[7-9],即

fi(x)=fi(x,y)=fi1(x)fi2(y)

(7)

常用的自相關(guān)函數(shù)型式有指數(shù)型(exp)與指數(shù)平方型(exp2)[1,7],其表達(dá)式分別為:

(8)

(9)

其中,Lh、Lv分別為水平、豎直方向上隨機(jī)場的相關(guān)長度(自相關(guān)距離);x1=(x1,y1),x2=(x2,y2),表示隨機(jī)場的空間位置。

1.2 KL展開法的誤差

理論上,(1)式中的M應(yīng)該取無窮大,但在保證計(jì)算精度的前提下,通常只截取前M項(xiàng)來提高計(jì)算效率。因此,如何定義隨機(jī)場的離散誤差,對于M正確選取及隨機(jī)場離散精度具有重要影響。

1.2.1 比率法

文獻(xiàn)[4-6]建議采用隨機(jī)場期望能的比率因子確定截?cái)囗?xiàng)數(shù),認(rèn)為比率因子越接近1時(shí)離散誤差越小,對應(yīng)的離散誤差(本文記為ε1)為:

(10)

其中,S為離散區(qū)域的面積。

1.2.2 逐點(diǎn)均方誤差法

當(dāng)采用KL展開法進(jìn)行隨機(jī)場的離散時(shí),逐點(diǎn)均方誤差[7]可表示為:

(11)

展開(11)式,可得:

(12)

由(12)式可得整個(gè)區(qū)域上逐點(diǎn)誤差的平均值(本文記為ε2)為:

(13)

其中,Ne為離散點(diǎn)的個(gè)數(shù)。

理論上,由(11)式可得逐點(diǎn)均方誤差ε(x)應(yīng)為0或正值。但計(jì)算表明,部分離散點(diǎn)處ε(x)可能為負(fù),這使得采用(13)式求解ε(x)的平均值時(shí),正、負(fù)誤差中和可能使求得的平均誤差ε2偏小,不能反映整個(gè)離散區(qū)域內(nèi)誤差的真實(shí)情況。因此,本文先對每點(diǎn)的離散誤差取絕對值,再對所有離散點(diǎn)誤差的絕對值取平均(本文記為ε3)來表示整個(gè)區(qū)域離散誤差的平均值,即

(14)

1.2.3 3種誤差的對比分析

上述計(jì)算誤差的3種方法中,ε(x)是逐點(diǎn)誤差,與離散點(diǎn)的具體位置有關(guān);ε1、ε2、ε3是離散區(qū)域上的平均誤差,它們能更好地代表整個(gè)離散區(qū)域上隨機(jī)場的離散精度,是決定M的重要依據(jù)。

ε1與ε2的關(guān)系為:當(dāng)隨機(jī)場離散的網(wǎng)格密度較大時(shí),(13)式中的求和符號將變?yōu)榉e分符號;將(12)式代入(13)式,由特征函數(shù)的正交性可知,ε1與ε2的計(jì)算公式相同,因此,當(dāng)隨機(jī)場離散的網(wǎng)格密度足夠大時(shí),ε2趨近于ε1,即ε1是ε2的極限形式。

ε2與ε3的關(guān)系為:若ε(x)全為非負(fù)數(shù)值,則ε2=ε3;否則,ε2<ε3。

2 算例分析

2.1 算例簡介

本文用2個(gè)算例來進(jìn)行計(jì)算分析。

例1 對一個(gè)桿件進(jìn)行隨機(jī)場離散,屬于一維空間問題。該桿件長l1=10 m。

例2 對一個(gè)矩形域進(jìn)行隨機(jī)場離散,屬于二維空間問題。該矩形區(qū)域的長度l1=14 m,高度l2=6 m。該離散區(qū)域?qū)?yīng)于平面應(yīng)變問題的地基承載力分析,其長度及高度方向分別表示土體中的水平及豎直方向。

2個(gè)算例對應(yīng)的各參數(shù)基準(zhǔn)值見表1所列。其中,對于二維空間問題,根據(jù)文獻(xiàn)[10-13]的研究結(jié)果,令長度、高度方向(即水平、豎直方向)的自相關(guān)距離各不相同,分別為30、3 m。

表1 隨機(jī)場的統(tǒng)計(jì)特征及離散條件基準(zhǔn)值

2.2 隨機(jī)場的逐點(diǎn)離散誤差

由(10)～(14)式可知,ε1只與特征值有關(guān);ε2、ε3只與特征值、特征函數(shù)有關(guān)。因此,當(dāng)離散區(qū)域大小、網(wǎng)格密度、展開項(xiàng)數(shù)、自相關(guān)函數(shù)類型及自相關(guān)距離等條件確定時(shí),隨機(jī)場的逐點(diǎn)離散誤差與平均誤差是定值,不隨隨機(jī)場離散次數(shù)的變化而變化。在表1基準(zhǔn)值條件下,例1與例2隨機(jī)場的逐點(diǎn)離散誤差如圖1所示。

由圖1可知,無論是一維問題還是二維問題,隨機(jī)場的逐點(diǎn)離散誤差都是在區(qū)域邊界上相對較大。由(11)式可知,ε(x)應(yīng)為非負(fù)數(shù)值。從圖1可以看出,由(12)式求得的逐點(diǎn)均方誤差ε(x)有部分值為負(fù)。產(chǎn)生這種現(xiàn)象的原因是:

(1) 在KL展開法中,λi是自相關(guān)函數(shù)ρ(x1,x2)的前M個(gè)由大到小排列的特征值;但fi(x)只是與λi對應(yīng)的特征函數(shù),它并非由大到小排列,這使得在只截取無窮級數(shù)的前M項(xiàng)作為隨機(jī)場離散值的近似值時(shí),若前M項(xiàng)中含有的特征函數(shù)值fi(x)較大,而相應(yīng)的λi又不是很小,則逐點(diǎn)均方誤差ε(x)就有可能為負(fù)。

(2) 由(2)式可知,特征函數(shù)值與正交基展開函數(shù)h(x)有關(guān),而后者又與標(biāo)準(zhǔn)Legendre多項(xiàng)式函數(shù)及離散點(diǎn)的具體位置x有關(guān)。例如,M=4,k為1～M時(shí),標(biāo)準(zhǔn)Legendre多項(xiàng)式Lk-1(u)=Lk-1((x-T)/a)的函數(shù)圖形如圖2所示。

從圖2可以看出,對于離散域邊界上的點(diǎn),Lk-1(u)恒為1;對于離散域邊界附近的點(diǎn),Lk-1(u)的絕對值也較大。因此,離散域邊界上及邊界附近,hk(x)的絕對值隨著k的增加而增加(見(5)式),導(dǎo)致特征函數(shù)fi(x)也較大(見(2)式),這使得離散區(qū)域邊界及其附近離散點(diǎn)上的逐點(diǎn)均方誤差ε(x)極易為負(fù)值。

圖1 隨機(jī)場的逐點(diǎn)離散誤差

圖2 k取不同值時(shí)標(biāo)準(zhǔn)Legendre多項(xiàng)式曲線

2.3 隨機(jī)場離散誤差的影響因素

以表1的參數(shù)值為基準(zhǔn)值,分別改變離散區(qū)域大小、離散網(wǎng)格密度、展開項(xiàng)數(shù)M、自相關(guān)函數(shù)類型及自相關(guān)距離,研究其對隨機(jī)場平均離散誤差ε1、ε2、ε3的影響。分析某個(gè)影響因素時(shí),只改變此影響因素的相應(yīng)值,而其他計(jì)算條件仍取表1的基準(zhǔn)值。

2.3.1 離散區(qū)域大小及離散網(wǎng)格密度的影響

分別取離散區(qū)域的邊長為基準(zhǔn)條件的1～3倍(對于二維算例,矩形區(qū)域的2條邊長等比例增加),再對2種離散區(qū)域,分別取隨機(jī)場網(wǎng)格邊長為0.25、0.50、1.00 m。離散區(qū)域大小、離散網(wǎng)格密度對隨機(jī)場離散誤差的影響分別如圖3、圖4所示。

圖3 隨機(jī)場離散誤差與離散區(qū)域大小的關(guān)系

圖4 隨機(jī)場離散誤差與網(wǎng)格密度的關(guān)系

圖3中,橫坐標(biāo)l1/Lh表示離散區(qū)域的相對大小。從圖3可以看出,無論是一維算例還是二維算例,3種誤差均隨著離散區(qū)域的增大而快速增加;ε1與ε2較為接近(ε1略大于ε2),ε3較大。圖4中,橫坐標(biāo)所示的網(wǎng)格邊長越小,則網(wǎng)格密度越大。從圖4可以看出,ε1與網(wǎng)格密度無關(guān);ε2隨網(wǎng)格密度的增加而有所增加,網(wǎng)格越密,ε2越接近于ε1;ε3隨網(wǎng)格密度的增加無明顯變化規(guī)律。圖4結(jié)果與1.2.3節(jié)的理論分析結(jié)果一致。

對比圖3、圖4的縱坐標(biāo)可以發(fā)現(xiàn),圖3縱坐標(biāo)的數(shù)值變化很大,圖4縱坐標(biāo)的數(shù)值變化很小。因此,離散區(qū)域大小對隨機(jī)場離散誤差的影響很大,網(wǎng)格密度對隨機(jī)場離散誤差的影響則很小,相比之下可以忽略不計(jì)。實(shí)際上,隨機(jī)場離散的KL展開法是一種級數(shù)展開法,其離散結(jié)果是連續(xù)的,與隨機(jī)場的網(wǎng)格密度無關(guān)。上述計(jì)算結(jié)果中,ε2、ε3與網(wǎng)格密度有關(guān),這是由于ε2、ε3只計(jì)算了特定離散點(diǎn)處的平均離散誤差(見(13)式、(14)式)。

2.3.2 展開項(xiàng)數(shù)的影響

在采用KL展開法進(jìn)行隨機(jī)場的離散時(shí),展開項(xiàng)數(shù)M是關(guān)系到離散精度的重要因素。M越大,離散誤差越小;但是,隨著M增加,隨機(jī)場離散的計(jì)算工作量也大幅增加,計(jì)算速度變慢。因此,需要確定合理的M?；鶞?zhǔn)值條件下,2個(gè)算例的隨機(jī)場離散誤差與M的關(guān)系如圖5所示。

圖5 隨機(jī)場離散誤差與M的關(guān)系

(1) 當(dāng)M較小時(shí),ε1、ε2明顯小于ε3,因此,根據(jù)ε1、ε2、ε3是否小于某個(gè)誤差限值來確定M時(shí),3種方法得到的M值并不相同。例如,設(shè)誤差限值為ε0=5%,對于例1,由ε1<ε0及ε2<ε0可得M=3,由ε3<ε0可得M=4;對于例2,由ε1<ε0可得M=4,由ε2<ε0可得M=3,由ε3<ε0可得M=4。

(2)ε1、ε2、ε3總體上隨M增加而快速減小,但二維算例的ε2在M=3時(shí)明顯偏小,這是由于此條件下逐點(diǎn)誤差的過度正負(fù)中和而導(dǎo)致的,因此,ε2不適用于判斷M。由于ε1是ε2的極限形式,因此,ε1也不適用于判斷M。ε3能綜合反映離散域上各點(diǎn)誤差的算術(shù)平均值,是判斷M的較好判別標(biāo)準(zhǔn)。

2.3.3 自相關(guān)函數(shù)類型、自相關(guān)距離的影響

自相關(guān)函數(shù)類型、自相關(guān)距離大小對隨機(jī)場離散誤差的影響如圖6所示。

圖6 隨機(jī)場離散誤差與自相關(guān)函數(shù)及自相關(guān)距離的關(guān)系

因?yàn)棣?不適用于判斷M,而ε1是ε2的極限形式,所以只繪了ε3與自相關(guān)函數(shù)類型、自相關(guān)距離大小的關(guān)系。由圖6可知,exp2對應(yīng)的離散誤差明顯小于exp對應(yīng)的誤差,因此對于一個(gè)具體的研究對象,如果對隨機(jī)場的自相關(guān)函數(shù)類型沒有明確的試驗(yàn)數(shù)據(jù),可以采用exp2,其對應(yīng)的離散誤差偏小;自相關(guān)函數(shù)類型相同時(shí),自相關(guān)距離越大,隨機(jī)場離散誤差越小,這是由于自相關(guān)距離越大時(shí),同樣條件下對應(yīng)的自相關(guān)函數(shù)值越大(見(8)式、(9)式),由(4)式得到的矩陣A的元素值越大,其對應(yīng)的特征值也越大,因此離散面積不變時(shí),隨機(jī)場的離散誤差就越小。

3 結(jié) 論

(1) 對于隨機(jī)場離散的KL展開法,可以采用3種方法求解離散域上的平均誤差。誤差ε1是誤差ε2的極限形式;誤差ε3可以合理反映逐點(diǎn)均方誤差為負(fù)的情況,是判斷隨機(jī)場展開項(xiàng)數(shù)M的一種較好的誤差判別標(biāo)準(zhǔn)。

(2) 離散區(qū)域大小增加時(shí),離散誤差ε1、ε2、ε3增加;離散網(wǎng)格密度增加時(shí),ε1不變,ε2趨近于ε1。

(3) 指數(shù)平方型自相關(guān)函數(shù)對應(yīng)的離散誤差明顯小于指數(shù)型自相關(guān)函數(shù)對應(yīng)的誤差。

(4) 自相關(guān)函數(shù)類型相同時(shí),自相關(guān)距離越大,隨機(jī)場離散誤差越小。