基于比例邊界有限元法計算應力強度因子的不確定量化分析

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

(1.防災減災湖北省重點實驗室(三峽大學),湖北 宜昌 443002;2.三峽大學 土木與建筑學院,湖北 宜昌 443002)

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

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

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

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

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

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

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

1.1 比例邊界有限元法

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

圖1 比例邊界有限元的坐標系統

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

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

(1)

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

應變表示為

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

(2)

式中,[B1(η)]和[B2(η)]為節點應變與位移的關系。

應力表示為

(3)

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

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

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

(4)

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

(5)

(6)

(7)

(8)

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

1.2 計算應力強度因子

求解式(4)后,可以獲得節點位移。通過位移場可以提取應力強度因子。位移場的漸進解表示為

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

(9)

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

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

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

O{ξ2}

(10)

位移的導數求解為

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

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

O{ξ}

(11)

其中,積分常數{c}表示為

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

(12)

在一個線單元的指定局部坐標ξ處,應力場表示為

(13)

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

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

(14)

奇異點周圍應力的漸近解表明,動力學中奇異應力和T-應力的定義與靜力學相同。圖2給出了基于SBFEM多邊形裂紋模型,其中笛卡爾坐標系和極坐標系的原點位于裂紋尖端。對于均質板中的裂紋和各向同性雙材料板中的界面裂紋,應力強度因子統一定義為

(15)

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

并且

(16)

式中:L為特征長度,應力強度因子的幅值與L無關;ε為兩種材料的振蕩系數,只與材料本身有關。

(17)

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

(18)

使用式(17),將奇異應力式(13)從局部坐標轉換為極坐標形式

(19)

根據式(15)中應力強度因子的形式,廣義應力強度因子表示為

(20)

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

1.3 蒙特卡羅模擬

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

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

(21)

對于靜力學問題略去后兩項的慣性力和阻尼力,邊界上的約束條件為

(22)

(23)

式中,{ε}為應變。則應力為

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

(24)

使用伽遼金方法,運用分部積分和邊界條件,可以得到帶有隨機參數的弱形式

(25)

使用比例邊界有限元法進行離散,得到帶有隨機參數的代數方程

(26)

式中:{d}為節點位移向量;{F}為荷載向量;[M]、[C]和[K]分別為全局質量矩陣、阻尼矩陣和剛度矩陣。值得注意的是,這里假設彈性模量作為隨機參數存在于每個單元剛度矩陣中,用隨機場f(x)來描述彈性模量,對于每個單元有

(27)

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

(28)

(29)

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

2 奇異值分解

根據初始樣本的波動范圍,引入m個隨機變量并計算其響應,將這些樣本點的響應[d(αm)]按列排列,構造快照矩陣[Ω][28]

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

(30)

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

(31)

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

(32)

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

(33)

式中:φ為基函數;χ為其權系數。本文選取高斯核函數作為基函數,如下

(34)

(35)

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

(36)

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

圖3 算法的流程圖

3 數值算例

3.1 帶斜裂紋的平面板

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

(37)

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

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

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

(38)

表1 不同變異系數條件下本文算法結果的精度

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

圖5 不同變異系數的預測值與解析解的比較

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

方案1:使用SBFEM計算1 000個隨機變量下的KI直接進行蒙特卡羅模擬。該計算需要重復1 000次。

方案2:對表1中的原始快照矩陣進行奇異值分解,截斷[Σ]矩陣,將[Σ]矩陣維數選擇為行與列的最小值。再進行蒙特卡羅模擬。

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

圖6給出了KI在不同變異系數下的統計特征。從圖中可以看出方案1與方案2較為吻合,方案3具有一定的誤差。這是由于方案3截斷了過多的奇異值,導致降階系統“失真”。為了避免系統失真,保留降階優勢,下面的算例都采用方案2的截斷方案。

圖6 不同方案下KI的統計特征

3.2 平面板中心彎曲裂紋

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

圖7 中心彎曲裂紋平面板

受均布荷載作用,計算采用平面應變假定。該問題具有解析解,應力強度因子為

(39)

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

表2 SBFEM與SVD-RBF進行隨機分析的計算成本

圖8 應力強度因子的算法解與解析解的比較

3.3 帶孔洞的雙裂紋矩形板

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

(a) 幾何尺寸

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

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

表3 SBFEM與SVD-RBF進行隨機分析的計算成本

圖10 無量綱動應力強度因子KI的統計特征

圖11 無量綱動應力強度因子KII的統計特征

圖12 本文算法解與參考解的動應力強度因子對比

3.4 帶有裂紋的重力壩-地基系統的地震響應分析

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

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

采用多邊形單元進行離散,劃分93個多邊形單元458節點,其網格如圖14所示。其中截斷邊界采用11個三節點線單元離散,相似中心選取為壩基的中點處,采用基于連分式算法的高階透射邊界模擬壩基無限域。選取動態應力強度因子KI作為響應,計算時間為10 s,步長Δt=0.01 s。

圖14 SBFEM網格

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

圖15 歸一化動應力強度因子KI的統計特征

圖16 Er=12 GPa時本文算法解與參考解的結果對比

4 結 論

本文構建了一種新的使用比例邊界有限元法計算應力強度因子的蒙特卡羅模擬框架,并得到如下結論:

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

(2) 采用蒙特卡羅模擬計算響應的統計特征,量化不確定參數對結構響應的影響,并使用模型降階法對傳統的蒙特卡羅模擬進行改進,直接降低物理問題的自由度,快速獲得新的結構響應,并討論了奇異值矩陣的合理截斷。同時,考慮了不同變異系數對結構響應的影響,對大的變異系數(0.4)的結果也有較高的精度。

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

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