劉雪瑩, 譚曉慧, 費鎖柱, 侯曉亮

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

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

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

1 隨機場離散的KL展開法

1.1 KL展開法的基本原理

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

(1)

其中,ξi為M個互為獨立的標準正態(tài)變量;M為隨機場的展開項數(shù)(截斷項數(shù));λi為自相關函數(shù)ρ(x1,x2)前M個由大到小排列的特征值;fi(x)為與λi對應的特征函數(shù)。

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

(2)

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

AD=DΛ

(3)

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

(4)

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

k=1,2,…,∞

(5)

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

(6)

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

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

(7)

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

(8)

(9)

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

1.2 KL展開法的誤差

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

1.2.1 比率法

文獻[4-6]建議采用隨機場期望能的比率因子確定截斷項數(shù),認為比率因子越接近1時離散誤差越小,對應的離散誤差(本文記為ε1)為:

(10)

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

1.2.2 逐點均方誤差法

當采用KL展開法進行隨機場的離散時,逐點均方誤差[7]可表示為:

(11)

展開(11)式,可得:

(12)

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

(13)

其中,Ne為離散點的個數(shù)。

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

(14)

1.2.3 3種誤差的對比分析

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

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

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

2 算例分析

2.1 算例簡介

本文用2個算例來進行計算分析。

例1 對一個桿件進行隨機場離散,屬于一維空間問題。該桿件長l1=10 m。

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

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

表1 隨機場的統(tǒng)計特征及離散條件基準值

2.2 隨機場的逐點離散誤差

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

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

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

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

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

圖1 隨機場的逐點離散誤差

圖2 k取不同值時標準Legendre多項式曲線

2.3 隨機場離散誤差的影響因素

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

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

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

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

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

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

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

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

在采用KL展開法進行隨機場的離散時,展開項數(shù)M是關系到離散精度的重要因素。M越大,離散誤差越小;但是,隨著M增加,隨機場離散的計算工作量也大幅增加,計算速度變慢。因此,需要確定合理的M?；鶞手禇l件下,2個算例的隨機場離散誤差與M的關系如圖5所示。

圖5 隨機場離散誤差與M的關系

(1) 當M較小時,ε1、ε2明顯小于ε3,因此,根據(jù)ε1、ε2、ε3是否小于某個誤差限值來確定M時,3種方法得到的M值并不相同。例如,設誤差限值為ε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時明顯偏小,這是由于此條件下逐點誤差的過度正負中和而導致的,因此,ε2不適用于判斷M。由于ε1是ε2的極限形式,因此,ε1也不適用于判斷M。ε3能綜合反映離散域上各點誤差的算術平均值,是判斷M的較好判別標準。

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

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

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

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

3 結 論

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

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

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

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