胡昊文,陳燈紅,王乾峰,胡記磊,駱 歡

(1.防災(zāi)減災(zāi)湖北省重點(diǎn)實(shí)驗(yàn)室(三峽大學(xué)),湖北 宜昌 443002;2.三峽大學(xué) 土木與建筑學(xué)院,湖北 宜昌 443002)

裂紋存在于眾多結(jié)構(gòu)中,例如混凝土結(jié)構(gòu)通常是帶裂縫工作的,巖塊表面存在著許多裂隙,等。斷裂問題一直都是計(jì)算固體力學(xué)的難點(diǎn)和重點(diǎn)。其中,應(yīng)力強(qiáng)度因子是預(yù)測載荷作用下結(jié)構(gòu)中裂紋的產(chǎn)生和擴(kuò)展的重要參數(shù)。分析斷裂問題時(shí),有限元法(finite element method,FEM)以其成熟的商用性和對不同邊界、不同荷載狀況下的普遍適用性成為求解斷裂問題最主要的數(shù)值方法。但傳統(tǒng)有限元法在處理斷裂問題時(shí)會(huì)面臨一些挑戰(zhàn)。比如三角形或四邊形單元不能準(zhǔn)確地計(jì)算裂紋尖端的奇異應(yīng)力,并且裂紋尖端的網(wǎng)格需要幾個(gè)數(shù)量級的加密。一些學(xué)者使用邊界元法[1-2](boundary element method,BEM)來處理斷裂問題。邊界元法只需離散結(jié)構(gòu)表面,將問題的維數(shù)減少一維。其作為一種半解析法,邊界元法能夠更準(zhǔn)確地計(jì)算應(yīng)力強(qiáng)度因子。對于復(fù)雜問題,邊界元法所需的基本解很難解析獲得甚至是不存在,同時(shí)積分方程還會(huì)產(chǎn)生奇異積分。這些困難限制了邊界元法的發(fā)展。

Song等[3-4]提出的比例邊界有限元法(scaled bountary finite element method,SBFEM),是一種新的半解析數(shù)值方法,它克服了傳統(tǒng)有限元法的一些局限性,并借鑒邊界元法的思想,只在結(jié)構(gòu)表面進(jìn)行離散,將問題維數(shù)減少一維,輻射邊界條件自然滿足,且不需要基本解。這些優(yōu)勢使其成為學(xué)者們所青睞的一種數(shù)值方法。該方法已成功應(yīng)用于無限域問題[5-8]、斷裂力學(xué)[9-12]、波傳播[13-14]、接觸問題[15]等多個(gè)領(lǐng)域,關(guān)于SBFEM最新研究詳見相關(guān)綜述[16]。當(dāng)處理線性斷裂力學(xué)時(shí),任何類型的奇異應(yīng)力都能解析表示,且不需要局部細(xì)化和富集函數(shù)[17]。此外,比例邊界有限元法多邊形單元的形函數(shù)包含奇異項(xiàng),可以準(zhǔn)確地表征應(yīng)力集中系數(shù)[18]、應(yīng)力強(qiáng)度因子[19]和T-應(yīng)力[20]。

目前成熟的數(shù)值算法,大都是基于確定性的分析。然而,在工程問題中存在諸多不確定性因素,例如結(jié)構(gòu)的幾何尺寸,材料自身的力學(xué)參數(shù),荷載的隨機(jī)性等等。若忽略這些因素的影響,可能會(huì)導(dǎo)致結(jié)構(gòu)失效,甚至引發(fā)災(zāi)難性的后果。綜合考慮不確定性因素,使數(shù)值分析更符合實(shí)際工程,不確定性量化[21](uncertainty quantification,UQ)分析成為近些年的熱點(diǎn)方向。

不確定量化本質(zhì)上是研究輸入的不確定造成的輸出不確定。針對這一特點(diǎn),學(xué)者們處理不確定量化分析的步驟主要分為三步。第一步:構(gòu)建物理模型。不確定量化的結(jié)果取決于高精度的數(shù)值模型。這一步通常采用成熟的確定性算法,諸如有限元法、邊界元法、比例邊界有限元法等。不確定量化與這些數(shù)值算法結(jié)合,產(chǎn)生了隨機(jī)有限元法、隨機(jī)邊界元法等經(jīng)典研究方法。第二步:描述不確定性參數(shù)。通常使用隨機(jī)變量或隨機(jī)場。例如攝動(dòng)法,將一個(gè)隨機(jī)函數(shù)在其均值附近進(jìn)行Taylor展開,然后對其進(jìn)行合理的截?cái)?多項(xiàng)式混沌展開,在隨機(jī)參數(shù)空間對精確解作一個(gè)有限階多項(xiàng)式展開,求解相關(guān)系數(shù)以獲得精確解的全部信息;蒙特卡羅模擬(Monte Carlo simulation,MCS)[22-23],假設(shè)隨機(jī)變量所遵循的概率分布,通過計(jì)算隨機(jī)抽樣的樣本響應(yīng),獲得統(tǒng)計(jì)特征(期望、標(biāo)準(zhǔn)差等)。第三步:將隨機(jī)變量或隨機(jī)場導(dǎo)入數(shù)值模型中,進(jìn)行不確定性傳播,得到高階矩、概率密度函數(shù)、累積密度函數(shù)等統(tǒng)計(jì)特征,用于評估不確定性參數(shù)對結(jié)構(gòu)響應(yīng)的影響。

描述不確定性參數(shù)步驟中,蒙特卡羅模擬方法以其思想簡單和廣泛適用成為很多學(xué)者首選的方法。此外,蒙特卡羅模擬的收斂速度只與隨機(jī)采樣次數(shù)N有關(guān),而與求解問題的維度無關(guān)。所以,在很多高維問題的求解過程中,蒙特卡羅方法往往是最合適的甚至是唯一的選擇。處理工程問題時(shí),單點(diǎn)響應(yīng)計(jì)算耗時(shí)嚴(yán)重,執(zhí)行大量采樣來確保蒙特卡羅模擬收斂的成本是無法負(fù)擔(dān)的,直接的蒙特卡羅模擬不得不與其他方法相結(jié)合,以節(jié)約計(jì)算成本,提高計(jì)算效率。Gu[24]從硬件方面減少硬件資源消耗、提高硬件計(jì)算的并行性和加速性能,從而減少蒙特卡羅模擬硬件方面的成本;Fishman[25]優(yōu)化蒙特卡羅模擬采樣策略,一定程度上減少了取樣個(gè)數(shù),降低了計(jì)算量;Ding等[26]使用模型降階算法加速蒙特卡羅模擬過程,分析結(jié)構(gòu)的高維不確定性。

本文構(gòu)建一種新穎的基于比例邊界有限元法求解斷裂問題的不確定量化分析方法。首先,采用比例邊界有限元法計(jì)算應(yīng)力強(qiáng)度因子(stress intensity factor,SIF);然后,使用蒙特卡羅模擬進(jìn)行不確定量化分析,并對直接的蒙特卡羅模擬進(jìn)行改進(jìn);最后采用奇異值分解,構(gòu)造低階的子空間,用徑向基函數(shù)對子空間進(jìn)行近似,通過子空間線性組合獲得新的結(jié)構(gòu)響應(yīng),實(shí)現(xiàn)基于MCS的快速不確定量化分析。

1 基于比例邊界有限元法的蒙特卡羅模擬

1.1 比例邊界有限元法

Wolf等詳細(xì)介紹了比例邊界有限元的推導(dǎo)過程,本文只給出一些必要的公式。在二維域中給出SBFEM所需的局部坐標(biāo),如圖1所示。選擇一個(gè)點(diǎn)O作為比例中心,對比例中心的唯一要求是S的整個(gè)邊界可視。徑向局部坐標(biāo)ξ,也被稱作比例因子,從比例中心指向邊界上任意一點(diǎn)。η為環(huán)向局部坐標(biāo)。這種由比例系數(shù)和單元局部坐標(biāo)定義的轉(zhuǎn)換關(guān)系稱為比例邊界轉(zhuǎn)化,并且這種轉(zhuǎn)化關(guān)系是唯一的。對于有限域,0≤ξ≤1,而對于無限域,1≤ξ≤∞。在(η,ξ)點(diǎn)處的位移可以用分段插值來近似

圖1 比例邊界有限元的坐標(biāo)系統(tǒng)

{u(ξ,η)}=[N(η)]{u(ξ)}=[N1(η)[I],

N2(η)[I],…]{u(ξ)}

(1)

式中:[N(η)]為在環(huán)向上的形函數(shù);u(ξ)為沿徑向線的位移函數(shù),僅與ξ有關(guān);[I]為一個(gè)2×2的單位矩陣。

應(yīng)變表示為

{ε(ξ,η)}=[B1(η)]{u(ξ)},ξ+

(2)

式中,[B1(η)]和[B2(η)]為節(jié)點(diǎn)應(yīng)變與位移的關(guān)系。

應(yīng)力表示為

(3)

式中,[D]是彈性矩陣。

在比例邊界坐標(biāo)中表示控制微分方程后,可以在徑向ξ上運(yùn)用伽遼金法或虛功原理,得到基于位移的二維比例邊界有限元法的基本方程

[E0]ξ2{u(ξ)},ξξ+([E0]-[E1]+[E1]T)ξ{u(ξ)},ξ-

(4)

式中,當(dāng)頻率ω=0,式(4)為靜力平衡方程。矩陣[E0]、[E1]、[E2]和[M0]是組裝邊界(ξ= 1)上計(jì)算的單元系數(shù)矩陣獲得,其形式為

(5)

(6)

(7)

(8)

值得說明的是,[E0]、[E1]和[E2]與常規(guī)有限元中靜力剛度陣相似,而[M0]類似于質(zhì)量陣。其中,[E0]、[M0]是正定矩陣,[E2]是半正定矩陣,而[E1]是非對稱。式(4)是Euler-Cauchy方程,可解析求解。

1.2 計(jì)算應(yīng)力強(qiáng)度因子

求解式(4)后,可以獲得節(jié)點(diǎn)位移。通過位移場可以提取應(yīng)力強(qiáng)度因子。位移場的漸進(jìn)解表示為

{u(ξ)}=[V11]ξ-[S11]{c}+O{ξ2}

(9)

式中:{c}為邊界積分常數(shù);O{ξn}為一個(gè)高階項(xiàng),在相同的ξn下趨于零,它包含高頻下慣性力的影響;[S]和[V]為式(4)轉(zhuǎn)化為一階常微分方程求解過程中Hamilton矩陣進(jìn)行Schur分解時(shí)產(chǎn)生的特征值和特征向量[27];[S11]為包含所有非正特征值;[V11]為對應(yīng)的特征值向量。代入到式(9)得

{u(ξ)}=[ψ(r)]{c(r)}+[ψ(s)]ξ-[S(s)]{c(s)}+

[ψ(T)]ξ{c(T)}+[ψn1]ξ-[Sn1]{cn1}+

O{ξ2}

(10)

位移的導(dǎo)數(shù)求解為

{u(ξ)},ξ=-[ψ(s)][S(s)]ξ-[S[s]]-[I]+

[ψ(T)]{c(T)}-[ψn1][Sn1]ξ-[Sn1]-[I]{cn1}+

O{ξ}

(11)

其中,積分常數(shù){c}表示為

{c}={{cn1};{c(T)};{c(s)};{c(r)}}

(12)

在一個(gè)線單元的指定局部坐標(biāo)ξ處,應(yīng)力場表示為

(13)

[B2(η)][ψ(s)]),

[B2(η)][ψ(T)])

(14)

奇異點(diǎn)周圍應(yīng)力的漸近解表明,動(dòng)力學(xué)中奇異應(yīng)力和T-應(yīng)力的定義與靜力學(xué)相同。圖2給出了基于SBFEM多邊形裂紋模型,其中笛卡爾坐標(biāo)系和極坐標(biāo)系的原點(diǎn)位于裂紋尖端。對于均質(zhì)板中的裂紋和各向同性雙材料板中的界面裂紋,應(yīng)力強(qiáng)度因子統(tǒng)一定義為

(15)

圖2 裂紋區(qū)域的離散(比例中心位于裂紋尖端)

并且

(16)

式中:L為特征長度,應(yīng)力強(qiáng)度因子的幅值與L無關(guān);ε為兩種材料的振蕩系數(shù),只與材料本身有關(guān)。

(17)

式中,r(θ)為從比例中心到沿角度(θ)徑向線的邊界的距離(圖2)。將ξ的矩陣冪函數(shù)改寫為極坐標(biāo)形式

(18)

使用式(17),將奇異應(yīng)力式(13)從局部坐標(biāo)轉(zhuǎn)換為極坐標(biāo)形式

(19)

根據(jù)式(15)中應(yīng)力強(qiáng)度因子的形式,廣義應(yīng)力強(qiáng)度因子表示為

(20)

通過式(20)可以很容易地計(jì)算任意(θ)角度下應(yīng)力強(qiáng)度因子,并且對于確定裂紋擴(kuò)展是非常重要的。

1.3 蒙特卡羅模擬

進(jìn)行蒙特卡羅模擬之前,首先簡要介紹帶有不確定性參數(shù)的控制方程。設(shè)整個(gè)分析域?yàn)棣?邊界為Γ(Γ=Γu+Γt),假定彈性模量E為不確定性參數(shù),即

{σ(E)}ij,j+{f}i-ρ{u(E)}i,tt-μ{u(E)}i,t=0

(21)

對于靜力學(xué)問題略去后兩項(xiàng)的慣性力和阻尼力,邊界上的約束條件為

(22)

(23)

式中,{ε}為應(yīng)變。則應(yīng)力為

{σ(E)}ij=[D(E)]ijkl{ε(E)}kl

(24)

使用伽遼金方法,運(yùn)用分部積分和邊界條件,可以得到帶有隨機(jī)參數(shù)的弱形式

(25)

使用比例邊界有限元法進(jìn)行離散,得到帶有隨機(jī)參數(shù)的代數(shù)方程

(26)

式中:7tpxvz5為節(jié)點(diǎn)位移向量;{F}為荷載向量;[M]、[C]和[K]分別為全局質(zhì)量矩陣、阻尼矩陣和剛度矩陣。值得注意的是,這里假設(shè)彈性模量作為隨機(jī)參數(shù)存在于每個(gè)單元?jiǎng)偠染仃囍?用隨機(jī)場f(x)來描述彈性模量,對于每個(gè)單元有

(27)

式中:[B]為位移微分矩陣;J為雅可比矩陣。將式(25)中的[B]轉(zhuǎn)化為相應(yīng)的系數(shù)矩陣[E],得到帶有不確定參數(shù)的比例邊界有限元方程。(對于靜力學(xué)問題,略去質(zhì)量陣和阻尼陣。)隨后基于蒙特卡羅模擬的思想,將隨機(jī)變量依次代入系統(tǒng)中進(jìn)行求解,通過式(28)和(29)計(jì)算統(tǒng)計(jì)特征,以評估不確定參數(shù)對結(jié)構(gòu)響應(yīng)的影響,即

(28)

(29)

式中:N為樣本個(gè)數(shù);{d(j)}為第j個(gè)樣本點(diǎn)的結(jié)構(gòu)響應(yīng)。蒙特卡羅模擬的收斂速度為O(N-1/2),N越大,蒙特卡羅模擬收斂效果越好,計(jì)算成本也就越高。為了降低計(jì)算成本,提高計(jì)算效率。第2章將使用模型降階法,對直接的MCS進(jìn)行加速求解。

2 奇異值分解

根據(jù)初始樣本的波動(dòng)范圍,引入m個(gè)隨機(jī)變量并計(jì)算其響應(yīng),將這些樣本點(diǎn)的響應(yīng)[d(αm)]按列排列,構(gòu)造快照矩陣[Ω][28]

[Ω]=[d(α1),d(α2),…,d(αm)]=

(30)

式中:[Ω]∈Rn×m;n為單個(gè)輸入變量下的響應(yīng)個(gè)數(shù);m為變量個(gè)數(shù)。使用奇異值分解法(singular value decomposition,SVD)將矩陣[Ω]分解

(31)

式中:r=min(m,n);[U]∈Rn×n、[V]∈Rm×m為正交矩陣;uj和vj分別為[ΩΩT]和[ΩTΩ]的特征向量。[Σ]∈n×m是一個(gè)對角矩陣,對角元素σj按從大到小排列。定義[ψj]=uj,[ζj(αi)]=σjvij,式(31)可以改寫為

(32)

式中:[ψj]為正交基;[ζj(αi)]為相應(yīng)的基底。使用式(32),系統(tǒng)響應(yīng)可以通過簡化表達(dá)的[ψj]和[ζj(αi)]的線性組合來表示,且其規(guī)模比初始模型要小得多。為了實(shí)現(xiàn)對任意隨機(jī)輸入變量的系統(tǒng)響應(yīng)的連續(xù)近似,應(yīng)用徑向基函數(shù)(radial basis function,RBF)來對降階子空間中的基底進(jìn)行插值

(33)

式中:φ為基函數(shù);χ為其權(quán)系數(shù)。本文選取高斯核函數(shù)作為基函數(shù),如下

(34)

(35)

值得注意的是,根據(jù)樣本點(diǎn)不同,上述方程組是無條件非奇異的。將式(33)代入式(32),式(32)改寫為

(36)

可以看出任何系統(tǒng)響應(yīng)都可以通過使用這種線性組合來近似,無須通過SBFEM進(jìn)行完整地計(jì)算。(當(dāng)然,快照矩陣的形成仍然需要完整的計(jì)算)。圖3給出了本文算法的流程圖。與直接MCS相比,奇異值分解降低了系統(tǒng)的自由度,并提供了一個(gè)低階模型。RBF計(jì)算相關(guān)系數(shù),將子空間線性組合獲得新的響應(yīng)。SVD-RBF提高了計(jì)算效率,降低了計(jì)算成本,這些優(yōu)勢一定程度上擴(kuò)展了MCS的適用性。

圖3 算法的流程圖

3 數(shù)值算例

3.1 帶斜裂紋的平面板

圖4考慮一個(gè)在y方向,受拉力作用的帶有斜裂紋的平面板。其解析解[30]為

(37)

圖4 帶斜裂紋的平面板模型

式中,P為y方向施加的載荷。a是裂紋一半的長度;β是垂直方向與裂紋之間的夾角。彈性模量E=1 Pa,泊松比v=0.3,矩形板長度L=1 m。

表1給出了三種不同的變異系數(shù)下形成快照矩陣的樣本點(diǎn)數(shù)。并引入R2來評估所提出的方法的準(zhǔn)確性,其計(jì)算公式為

(38)

表1 不同變異系數(shù)條件下本文算法結(jié)果的精度

為了更直觀地展示算法的精度,圖5給出了不同變異系數(shù)下的算法解和解析解的比較。從圖中可以看到,變異系數(shù)控制著變量波動(dòng)的范圍(變異系數(shù)越大,數(shù)據(jù)波動(dòng)越大)。在各種變異系數(shù)下,本文算法與解析解吻合較好。結(jié)合表1,可以得出變異系數(shù)越大,保持算法精度所需的初始樣本點(diǎn)也就越多。初始樣本點(diǎn)的數(shù)目決定著計(jì)算效率,所以平衡計(jì)算效率與計(jì)算精度是值得注意的。

圖5 不同變異系數(shù)的預(yù)測值與解析解的比較

隨機(jī)分析之前,會(huì)進(jìn)行奇異值分解,此時(shí)產(chǎn)生了[Σ]矩陣。該矩陣只有主對角線上存在元素,并且從大到小排列。這些元素稱為奇異值,衰減迅速。奇異值越大包含的系統(tǒng)信息就越多,所以截?cái)郲Σ]矩陣是必要的,以節(jié)省儲(chǔ)存空間,提升計(jì)算效率。仍使用算例1的初始樣本,并采取三種方案進(jìn)行隨機(jī)性分析。

方案1:使用SBFEM計(jì)算1 000個(gè)隨機(jī)變量下的KI直接進(jìn)行蒙特卡羅模擬。該計(jì)算需要重復(fù)1 000次。

方案2:對表1中的原始快照矩陣進(jìn)行奇異值分解,截?cái)郲Σ]矩陣,將[Σ]矩陣維數(shù)選擇為行與列的最小值。再進(jìn)行蒙特卡羅模擬。

方案3:在方案2的基礎(chǔ)上,[Σ]矩陣維數(shù)選擇為σcat=10-1σmax。

圖6給出了KI在不同變異系數(shù)下的統(tǒng)計(jì)特征。從圖中可以看出方案1與方案2較為吻合,方案3具有一定的誤差。這是由于方案3截?cái)嗔诉^多的奇異值,導(dǎo)致降階系統(tǒng)“失真”。為了避免系統(tǒng)失真,保留降階優(yōu)勢,下面的算例都采用方案2的截?cái)喾桨浮?/p>

圖6 不同方案下KI的統(tǒng)計(jì)特征

3.2 平面板中心彎曲裂紋

本節(jié)考慮一個(gè)中心帶有彎曲裂紋的平面板問題[32]。其模型如圖7所示,半徑和中心裂紋角由r和θ定義。本例中,假設(shè)r/w=0.1,取r=4.25 m,上下邊界

圖7 中心彎曲裂紋平面板

受均布荷載作用,計(jì)算采用平面應(yīng)變假定。該問題具有解析解,應(yīng)力強(qiáng)度因子為

(39)

以應(yīng)力強(qiáng)度因子KI、KII和節(jié)點(diǎn)位移作為響應(yīng)構(gòu)造快照矩陣,并對應(yīng)力強(qiáng)度因子歸一化處理。選取初始裂紋角度作為隨機(jī)輸入變量,且服從期望μ= π/6,和標(biāo)準(zhǔn)差為σ=2的正態(tài)分布。設(shè)置分布參數(shù)后,使用隨機(jī)數(shù)生成器,隨機(jī)生成100個(gè)隨機(jī)角度變量,代入到解析解中進(jìn)行計(jì)算,作為蒙特卡羅模擬的參考解。圖8給出了本文算法與解析解的應(yīng)力強(qiáng)度因子KI,KII的對比。從圖中的結(jié)果來看,本文的算法與解析解的結(jié)果吻合得很好,表明算法具有較高的精度。值得注意的是,降階算法只采用了35個(gè)初始樣本點(diǎn)來形成快照矩陣,即只使用了35個(gè)樣本點(diǎn),獲得了100個(gè)該變量下的響應(yīng)結(jié)果,也說明了該算法在計(jì)算效率上具有優(yōu)勢。表2具體給出了該算法與基于SBFEM的蒙特卡羅模擬的統(tǒng)計(jì)特征和計(jì)算成本。表2中的快照矩陣由應(yīng)力強(qiáng)度因子和98個(gè)節(jié)點(diǎn)位移構(gòu)成,自由度的計(jì)算方式為變量個(gè)數(shù)乘以模型自由度。從表2可以看出,本文算法降低了系統(tǒng)自由度,加速了傳統(tǒng)的基于SBFEM的蒙特卡羅模擬過程,提高了計(jì)算效率。

表2 SBFEM與SVD-RBF進(jìn)行隨機(jī)分析的計(jì)算成本

圖8 應(yīng)力強(qiáng)度因子的算法解與解析解的比較

3.3 帶孔洞的雙裂紋矩形板

本節(jié)考慮一個(gè)帶有圓孔的矩形板,設(shè)圓孔半徑r為不確定變量,研究其對動(dòng)態(tài)應(yīng)力強(qiáng)度因子(dynamic stress intensity factor,DSIF)的影響。矩形板的寬度為2b=30 mm,長度為2h=60 mm,包含一個(gè)半徑為r的孔,在圓孔處有兩條與水平方向呈30°夾角的等長度裂紋,如圖9(a)所示。其材料特性分別為:密度ρ=5×10-6kg/mm3、泊松比v=0.3,剪切模量G=76.923 GPa。t=0時(shí),板的上下兩端受到均布的階躍荷載p(t)=PH(t)的作用,其中P是荷載的幅值。采用平面應(yīng)變假定。

(a) 幾何尺寸

使用多邊形比例邊界有限元法分析時(shí),將該模型劃分為8個(gè)多邊形單元,如圖9(b)所示。多邊形C1和C2是包含兩個(gè)裂紋的單元,并且將比例中心設(shè)置為裂紋的尖端。對于其他單元,都采用八節(jié)點(diǎn)的高階單元進(jìn)行離散,用正方形來表示每個(gè)單元兩端的節(jié)點(diǎn)。網(wǎng)格劃分完畢后,進(jìn)行隨機(jī)性分析。將圓孔半徑設(shè)為不確定變量,且滿足均值為μ=3.6,標(biāo)準(zhǔn)差為σ=0.2的正態(tài)分布。使用隨機(jī)數(shù)生成器產(chǎn)生50個(gè)滿足該分布的隨機(jī)變量,導(dǎo)入到模型中進(jìn)行計(jì)算,獲得50個(gè)變量下的動(dòng)態(tài)應(yīng)力強(qiáng)度因子作為蒙特卡羅模擬的參考解。其中,計(jì)算動(dòng)態(tài)應(yīng)力強(qiáng)度因子的總時(shí)間為20 μs,計(jì)算時(shí)間步長Δt=0.1 μs。

圖10和11給出了降階算法與蒙特卡羅模擬得出的統(tǒng)計(jì)特征的對比。從圖中的結(jié)果來看,本文的算法與蒙特卡羅模擬結(jié)果吻合得很好。值得注意的是,降階算法只采用了17個(gè)初始樣本點(diǎn)來形成快照矩陣。意味著降階算法只使用了17個(gè)樣本點(diǎn),獲得了50個(gè)蒙特卡羅模擬的結(jié)果。為了使數(shù)據(jù)更具說服性,圖12給出了r=3.75 mm時(shí),算法的結(jié)果與參考解的對比。其中,這里r=3.75 mm是算法預(yù)測值,而非初始樣本點(diǎn)。從圖12可以看出,算法預(yù)測的結(jié)果與參考解基本保持一致。此外,快照矩陣還使用奇異值分解進(jìn)行降維,進(jìn)一步節(jié)省儲(chǔ)存空間。為了更清晰地說明算法的優(yōu)越性,表3給出了基于SBFEM分析斷裂問題的框架下,本文算法與直接的蒙特卡羅模擬的計(jì)算成本的對比。總的來說,本文算法提高了傳統(tǒng)蒙特卡羅模擬的計(jì)算效率,拓展了蒙特卡羅模擬的適用性。

表3 SBFEM與SVD-RBF進(jìn)行隨機(jī)分析的計(jì)算成本

圖10 無量綱動(dòng)應(yīng)力強(qiáng)度因子KI的統(tǒng)計(jì)特征

圖11 無量綱動(dòng)應(yīng)力強(qiáng)度因子KII的統(tǒng)計(jì)特征

圖12 本文算法解與參考解的動(dòng)應(yīng)力強(qiáng)度因子對比

3.4 帶有裂紋的重力壩-地基系統(tǒng)的地震響應(yīng)分析

本節(jié)研究巖石地基對大壩動(dòng)力斷裂的影響。大壩的幾何尺寸如圖13所示。大壩材料為混凝土,其力學(xué)參數(shù)為:彈性模量Ec=31.0 GPa,泊松比vc=0.15,密度ρc=2 643 kg/m3;下部是巖石地基,部分力學(xué)參數(shù)vr=0.25,ρr=2 643 kg/m3。依照文獻(xiàn)[33],選取初始裂紋的高程為59.25 m,初始長度為1 m,輸入的地震動(dòng)為Koyna地震波。設(shè)定巖基的彈性模量為不確定性參數(shù),其波動(dòng)范圍為[10～15],步長為0.1選取50個(gè)Er作為蒙特卡羅模擬的結(jié)果。整個(gè)計(jì)算采用平面應(yīng)變假定。

圖13 Koyna大壩幾何模型 (m)

采用多邊形單元進(jìn)行離散,劃分93個(gè)多邊形單元458節(jié)點(diǎn),其網(wǎng)格如圖14所示。其中截?cái)噙吔绮捎?1個(gè)三節(jié)點(diǎn)線單元離散,相似中心選取為壩基的中點(diǎn)處,采用基于連分式算法的高階透射邊界模擬壩基無限域。選取動(dòng)態(tài)應(yīng)力強(qiáng)度因子KI作為響應(yīng),計(jì)算時(shí)間為10 s,步長Δt=0.01 s。

圖14 SBFEM網(wǎng)格

圖15給出了本文算法和直接蒙特卡羅模擬獲得的統(tǒng)計(jì)特征的對比圖。從圖15可以看到,兩種方法得到的統(tǒng)計(jì)特征基本一致。值得注意的是,算法只使用了11個(gè)初始樣本點(diǎn)形成快照矩陣,且采用方案2的策略對快照矩陣進(jìn)行降階處理。大大節(jié)省了儲(chǔ)存空間,提高了計(jì)算效率。為了進(jìn)一步說明算法的準(zhǔn)確性,圖16給出了Er=12 GPa算法的結(jié)果,與已有文獻(xiàn)結(jié)果[33]對比,這里的Er=12 GPa是算法預(yù)測得到結(jié)果而非初始樣本點(diǎn)。從圖16可以看出,算法解與參考解吻合得很好。KI的峰值發(fā)生在4 s處,在彈性模量Er的波動(dòng)范圍內(nèi),峰值處KI從8 641 527.9 Pa·m0.5波動(dòng)到7 023 890.08 Pa·m0.5相對波動(dòng)范圍為18%,可見巖基參數(shù)的不確定性對動(dòng)態(tài)應(yīng)力強(qiáng)度因子有較大影響,因此研究不確定參數(shù)對結(jié)構(gòu)響應(yīng)的影響是很有意義的。

圖15 歸一化動(dòng)應(yīng)力強(qiáng)度因子KI的統(tǒng)計(jì)特征

圖16 Er=12 GPa時(shí)本文算法解與參考解的結(jié)果對比

4 結(jié) 論

本文構(gòu)建了一種新的使用比例邊界有限元法計(jì)算應(yīng)力強(qiáng)度因子的蒙特卡羅模擬框架,并得到如下結(jié)論:

(1) 采用多邊形比例邊界有限元法分析斷裂問題時(shí),將比例中心設(shè)置在裂紋尖端,使裂紋尖端不需要網(wǎng)格細(xì)化,可直接提取應(yīng)力強(qiáng)度因子,為不確定量化提供了高精度的計(jì)算模型。

(2) 采用蒙特卡羅模擬計(jì)算響應(yīng)的統(tǒng)計(jì)特征,量化不確定參數(shù)對結(jié)構(gòu)響應(yīng)的影響,并使用模型降階法對傳統(tǒng)的蒙特卡羅模擬進(jìn)行改進(jìn),直接降低物理問題的自由度,快速獲得新的結(jié)構(gòu)響應(yīng),并討論了奇異值矩陣的合理截?cái)?。同時(shí),考慮了不同變異系數(shù)對結(jié)構(gòu)響應(yīng)的影響,對大的變異系數(shù)(0.4)的結(jié)果也有較高的精度。

(3) 數(shù)值算例表明,對于不同的變量和不同的物理模型,使用SVD-RBF加速M(fèi)CS具有良好的精度。

并且能顯著降低計(jì)算成本,提升了傳統(tǒng)蒙特卡羅模擬的適用性。