SBFEM 與非局部宏微觀損傷模型相結(jié)合的準(zhǔn)脆性材料開裂模擬1)

杜成斌 黃文倉 江守燕

(河海大學(xué)力學(xué)與材料學(xué)院,南京 211100)

引言

混凝土是一種典型的準(zhǔn)脆性材料,被廣泛應(yīng)用于土木和水利工程中.混凝土結(jié)構(gòu)的破壞通常是內(nèi)部微裂紋的萌生、發(fā)展貫通、形成宏觀裂紋并擴(kuò)展而引起結(jié)構(gòu)失效的過程[1].因此,研究混凝土結(jié)構(gòu)中裂紋的產(chǎn)生和擴(kuò)展的過程至關(guān)重要.自20 世紀(jì)20年代開始,科研工作者們對固體材料破壞機(jī)理進(jìn)行了大量深入研究,斷裂力學(xué)[2-4]與損傷力學(xué)[5-7]得到了快速發(fā)展.

采用傳統(tǒng)斷裂力學(xué)方法時需要在斷裂過程區(qū)劃分足夠細(xì)的網(wǎng)格,才能保證應(yīng)力強(qiáng)度因子計算的準(zhǔn)確性,且在模擬裂紋擴(kuò)展時需要復(fù)雜的網(wǎng)格重劃分[8],導(dǎo)致計算難度增大、效率降低.而采用傳統(tǒng)損傷力學(xué)進(jìn)行裂紋模擬時,僅需要在損傷過程區(qū)進(jìn)行網(wǎng)格加密,無需復(fù)雜的網(wǎng)格重劃分.損傷局部帶的模擬一直是計算力學(xué)研究的熱點.但傳統(tǒng)的局部損傷模型會引起應(yīng)變局部化,表現(xiàn)出明顯的網(wǎng)格敏感性[9].為了彌補(bǔ)局部模型的不足,非局部模型[9-11]逐漸發(fā)展起來.這些非局部模型使用非局部變量作為損傷驅(qū)動變量,一定程度上減輕了網(wǎng)格敏感性,但在非局部變量的正則化過程中容易產(chǎn)生錯誤的結(jié)果,導(dǎo)致不切實際的損傷啟動與擴(kuò)展[12-13].并且一些非局部公式在軟化階段會不可避免地出現(xiàn)應(yīng)力鎖死的問題[14],即在損傷區(qū)域兩側(cè)出現(xiàn)不該有的應(yīng)力傳遞.

近期,文獻(xiàn)[15-16]受近場動力學(xué)[17]基本思想與統(tǒng)一相場模型[18-19]的啟發(fā),提出了一種非局部宏微觀損傷模型.該模型引入物質(zhì)點對的概念,用物質(zhì)點對之間的物質(zhì)鍵變形來定義微細(xì)觀損傷,而宏觀尺度上的拓?fù)鋼p傷是由一定范圍內(nèi)微細(xì)觀損傷的累積引起的.通過引入能量退化函數(shù),將拓?fù)鋼p傷嵌入損傷力學(xué)的本構(gòu)關(guān)系中.該損傷模型不存在網(wǎng)格敏感性與應(yīng)力鎖死的問題[15],為準(zhǔn)脆性材料的開裂模擬提供了一種新思路.文獻(xiàn)[15]是在有限元框架下計算的,為了保證拓?fù)鋼p傷的計算精度,需要在損傷過程區(qū)劃分足夠精細(xì)的網(wǎng)格,這導(dǎo)致劃分網(wǎng)格與數(shù)值計算的時間成本大大增加.因此,采用更高效的網(wǎng)格劃分方案與數(shù)值計算方法是非常有必要的.

比例邊界有限元法(scaled boundary finite element method,SBFEM)是由Song 和Wolf[20]提出的一種新型半解析數(shù)值計算方法.該方法在環(huán)向上具有有限元精度的數(shù)值解,在徑向上具有解析解,在求解時需將求解域離散為一個或多個子域,且僅需采用一維線單元對每個子域的邊界進(jìn)行離散,使計算維度降低一維.將其與高效的四叉樹網(wǎng)格離散技術(shù)結(jié)合,不僅避免了懸掛節(jié)點的問題,還有效提高了計算效率[21].目前,基于四叉樹網(wǎng)格的比例邊界有限元已被成功應(yīng)用于動力分析[22]、裂紋擴(kuò)展[23]、彈塑性[24]、非局部損傷[25-26]等問題中.但總體而言,關(guān)于比例邊界有限元在損傷方面的研究還較少,有必要進(jìn)一步的探索研究.

本文將比例邊界有限元與非局部宏微觀損傷模型相結(jié)合,充分利用比例邊界有限元網(wǎng)格允許存在懸掛節(jié)點的優(yōu)勢,采用四叉樹網(wǎng)格離散技術(shù)實現(xiàn)快速、高質(zhì)量的多級網(wǎng)格劃分與尺寸過渡.將四叉樹網(wǎng)格六種單元的計算信息預(yù)先計算并存儲,計算時直接調(diào)用.以子域的比例中心作為物質(zhì)點,根據(jù)影響域范圍內(nèi)比例中心之間的相對變形,實現(xiàn)從物質(zhì)鍵微細(xì)觀損傷到宏觀拓?fù)鋼p傷的計算.通過兩個典型算例,說明本文方法可用于準(zhǔn)脆性材料的準(zhǔn)靜態(tài)開裂模擬,不存在網(wǎng)格敏感性問題,并有效提高了計算效率.

1 比例邊界有限元法基本理論

與常規(guī)有限元法類似,比例邊界有限元法在求解時需要將求解域離散為一個或多個子域,每個子域內(nèi)部需要設(shè)置一個比例中心,比例中心的選取必須滿足子域邊界上任意點都可見的條件[27].分別對每個子域求解后,即可得到整體的數(shù)值結(jié)果.由于該方法的相關(guān)文獻(xiàn)豐富,本文僅給出二維比例邊界有限元法的關(guān)鍵方程,關(guān)于更詳細(xì)的推導(dǎo)和解釋,可以參考文獻(xiàn)[21,27].

圖1 所示為一比例邊界有限元子域,比例中心O設(shè)置在子域的幾何中心處,子域邊界采用一維線單元離散.以O(shè)為原點建立局部坐標(biāo)系 (ξ,η),ξ 為徑向坐標(biāo),在比例中心O處 ξ=0,在邊界處 ξ=1,η 為環(huán)向坐標(biāo),其取值范圍為 [ -1,1] .對于子域內(nèi)任意一點的笛卡爾坐標(biāo) (x,y) 與局部坐標(biāo) (ξ,η) 的變換關(guān)系為

圖1 比例邊界有限元子域示意圖Fig.1 Subdomain of a scaled boundary finite element

式中 (x0,y0) 為比例中心坐標(biāo),xb=[x1,x2]T與yb=[y1,y2]T為邊界線單元的節(jié)點坐標(biāo),N(η) 為邊界線單元的插值形函數(shù).

子域內(nèi)任意一點 (ξ,η) 的位移可以表示為

式中,u(ξ) 為徑向位移函數(shù),可以通過求解式(3)的比例邊界有限元位移控制方程得到

該方程是關(guān)于 ξ 的二階齊次微分方程.E0,E1,E2為與材料參數(shù)和幾何形狀有關(guān)的系數(shù)矩陣,其形式為

式中,D為彈性矩陣,|J| 為 雅克比矩陣行列式,B1(η)和B2(η)為應(yīng)變位移轉(zhuǎn)換矩陣,具體形式見文獻(xiàn)[21].

式(3)為歐拉-柯西方程,可以解析求解.引入輔助變量X(ξ),可將方程轉(zhuǎn)換為一階齊次微分方程.令

則有

其中,q(ξ) 為與u(ξ) 對應(yīng)的節(jié)點力向量,Z為Hamilton矩陣,有

將Hamilton 矩陣進(jìn)行特征值分解,對于每個子域都滿足下式關(guān)系

式中,Λ=diag(λ1,λ2,···,λn)為Z矩陣正特征值組成的對角矩陣,Φu=[?u1,?u2,···,?un] 和Φq=[?q1,?q2,···,?qn]分別為位移模態(tài)和應(yīng)力模態(tài),n為子域的自由度數(shù).

通過模態(tài)疊加,式(3)所示二階齊次微分方程的解可以表示為

其中,c為與Z矩陣特征值相對應(yīng)的積分常數(shù),可以通過子域邊界節(jié)點位移ub=u(ξ=1) 求解得到

在比例中心處,只存在兩種剛體位移模態(tài)[21],而剛體位移模態(tài)由第n-1 階與第n階位移模態(tài)控制,故在比例中心的位移可以表示為

對于有界域,子域的剛度矩陣表示為[25]

對于彈性問題,將所有子域的剛度矩陣組裝形成整體剛度矩陣K,得到整體平衡方程

其中,U為整體位移向量,P為整體荷載向量.

2 非局部宏微觀損傷模型及其在SBFEM中的實現(xiàn)

非局部宏微觀損傷模型是一種基于幾何變形驅(qū)動的兩尺度損傷模型[28],該模型結(jié)合了近場動力學(xué)中物質(zhì)點對的基本思想與統(tǒng)一相場模型中能量退化函數(shù)的概念.為了便于理解,下面簡單介紹該模型的基本概念與思想,并將其嵌入比例邊界有限元的基本框架中,構(gòu)建基于比例邊界有限元框架的非局部宏微觀損傷模型.有關(guān)非局部宏微觀損傷模型更詳細(xì)的介紹可參考文獻(xiàn)[15-16].

2.1 微細(xì)觀損傷與拓?fù)鋼p傷

假設(shè)存在一連續(xù)體B,x與x′為B中的兩個物質(zhì)點.x與x′組成的物質(zhì)點對記為(x,x′),物質(zhì)點對(x,x′)之間建立的聯(lián)系稱為物質(zhì)鍵.根據(jù)圖2 所示的變形幾何,可定義物質(zhì)鍵的正伸長率為

圖2 物質(zhì)點對及其變形Fig.2 Deformation of material point pair

其中,ξ=x′-x為兩物質(zhì)點的初始空間差,u和u′分別為物質(zhì)點x和x′的位移,| |?|| 表 示向量長度,〈 ?〉 為麥克勞林算符,定義為 〈 α〉=max(0,α) .

當(dāng)物質(zhì)鍵的正伸長率超過材料的臨界伸長率λc時,物質(zhì)鍵開始發(fā)生不可逆的損傷,該過程與加載歷史參數(shù)有關(guān).為了表征物質(zhì)鍵的損傷程度,引入微細(xì)觀損傷函數(shù)ω∈[0,1],ω=0 表示無損狀態(tài),ω=1表示全損狀態(tài).顯然 ω 是加載歷史參數(shù)的函數(shù),可取為[15]

式中,γ 為材料參數(shù),與材料的微觀結(jié)構(gòu)有關(guān),可通過標(biāo)準(zhǔn)試件試驗確定.并且 γ 表征損傷發(fā)展的速度,γ的取值越大,損傷發(fā)展越快,脆性性質(zhì)越明顯,圖3給出了不同取值的 γ 對應(yīng)的微細(xì)觀損傷演化曲線.加載歷史參數(shù) κ 為物質(zhì)鍵的歷史最大超越伸長率,其定義為[15-16]

圖3 微細(xì)觀損傷演化曲線Fig.3 Microscopical damage evolution curve

式中,t為到達(dá)當(dāng)前荷載步的總時間,τ為t時刻之前某時刻所對應(yīng)荷載步的時間.通常取 λc=ft/E,ft為材料抗拉強(qiáng)度,E為彈性模量.

從宏觀上看,物質(zhì)點x的損傷狀態(tài)與其周圍一定范圍內(nèi)所連接的物質(zhì)鍵損傷有關(guān).宏觀損傷的發(fā)展代表著該點材料拓?fù)浒l(fā)生了改變[28],引入拓?fù)鋼p傷函數(shù) Ω (x) 來表征宏觀損傷程度.假設(shè)物質(zhì)點x的影響域半徑為l,影響域表示為Dl(x) .則物質(zhì)點x的拓?fù)鋼p傷可定義為影響域內(nèi)物質(zhì)鍵微細(xì)觀損傷的加權(quán)平均[15]

式中,影響域半徑l與材料的內(nèi)稟特征長度有關(guān)[9],而內(nèi)稟特征長度與材料的不均勻性相關(guān)(如骨料分布、骨料最大粒徑),可通過試驗來確定.φ (x,x′) 為權(quán)函數(shù),應(yīng)滿足非負(fù)性、定域性、歸一性[17].在本文中,取權(quán)函數(shù)為

式中,v ol[?] 表示體積測度.

后面算例表明,與常規(guī)的積分型非局部模型方法相比,本文的損傷模型計算可得到更準(zhǔn)確的局部開裂損傷帶,結(jié)果更為合理.

2.2 能量退化函數(shù)

物質(zhì)鍵發(fā)生損傷而引起拓?fù)鋼p傷發(fā)展的同時,還伴隨著能量耗散.引入表征物質(zhì)點儲能能力退化的能量退化因子[29]g,顯然g與拓?fù)鋼p傷的發(fā)展程度有關(guān),因此g是拓?fù)鋼p傷的函數(shù),記作能量退化函數(shù)g=g(Ω).g(Ω) 應(yīng)滿足g(0)=1,g(1)=0,由于損傷總是耗能的,則 dg(Ω)/dΩ ≤0 .對于經(jīng)典的連續(xù)介質(zhì)損傷力學(xué),g為 1 -Ω 的形式[5].應(yīng)當(dāng)指出能量退化函數(shù)g(Ω) 是區(qū)間 [ 0,1] 上的凸函數(shù),其定性的物理解釋可參考文獻(xiàn)[15].在統(tǒng)一相場模型[18]的能量退化函數(shù)的基礎(chǔ)上,采用兩參數(shù)的有理分式作為非局部宏微觀損傷模型中的能量退化函數(shù)[15],其形式為

式中,p和q均為退化參數(shù),其取值刻畫了材料的宏觀脆性程度,可以認(rèn)為是材料的物性參數(shù),可以通過標(biāo)準(zhǔn)試件試驗確定[16,30].

2.3 非局部宏微觀損傷模型在SBFEM 中的實現(xiàn)

如圖4 所示,在比例邊界有限元中,以子域的比例中心作為物質(zhì)點.首先根據(jù)式(11)計算比例中心的位移,從而通過式(14)計算比例中心之間所連物質(zhì)鍵的正伸長率,再根據(jù)式(15)對各物質(zhì)鍵的微細(xì)觀損傷進(jìn)行評估,最后通過式(17)對影響域內(nèi)物質(zhì)鍵的微細(xì)觀損傷進(jìn)行加權(quán)平均得到拓?fù)鋼p傷.如果在損傷過程區(qū)使用足夠小的子域,可假定損傷在一個子域內(nèi)是均勻分布的,即該子域的損傷程度可由該子域內(nèi)的物質(zhì)點所計算的拓?fù)鋼p傷值來表征.在經(jīng)典的連續(xù)介質(zhì)損傷力學(xué)中,本構(gòu)關(guān)系為[5]

其中,σ 為應(yīng)力張量,ε 為應(yīng)變張量.在非局部宏微觀損傷模型中,通過能量退化函數(shù)g=g(Ω) 的引入,建立起拓?fù)鋼p傷與應(yīng)力應(yīng)變之間的聯(lián)系,則本構(gòu)關(guān)系轉(zhuǎn)變?yōu)?/p>

式中,DD=gD為彈性損傷矩陣,將其代入式(4),可得到子域在損傷狀態(tài)下的系數(shù)矩陣

將式(22)代入式(7)的Z矩陣中,根據(jù)矩陣的性質(zhì),相應(yīng)的特征值、位移模態(tài)和應(yīng)力模態(tài)顯然與無損狀態(tài)下相同.再將式(22)代入式(12),得到子域在損傷狀態(tài)下的剛度矩陣

在損傷計算過程中,將所有的子域剛度矩陣進(jìn)行組裝,即可得到損傷狀態(tài)下的整體平衡方程

其中,KD為損傷狀態(tài)下的整體剛度矩陣.該方程為非線性方程,可以通過牛頓法、擬牛頓法、顯式迭代法、弧長法等方法求解.在本文中,由于下降段是裂紋模擬過程中的重要階段,故采用弧長法進(jìn)行求解,關(guān)于弧長法的具體概念與公式可參考文獻(xiàn)[31].

3 四叉樹網(wǎng)格

四叉樹是一種被廣泛使用的高效分層數(shù)據(jù)結(jié)構(gòu),特別適用于網(wǎng)格的快速劃分與網(wǎng)格尺寸快速平滑過渡.利用四叉樹離散技術(shù),可以高效、高質(zhì)量地建立起數(shù)值計算的多級網(wǎng)格模型.但在不同尺寸四叉樹網(wǎng)格之間存在懸掛節(jié)點,對于常規(guī)有限元無法直接使用四叉樹網(wǎng)格進(jìn)行計算,通常需要對網(wǎng)格進(jìn)行特殊處理.而SBFEM 可以將懸掛節(jié)點所在邊界視為兩個線單元,作為多邊形子域進(jìn)行計算,直接避免了懸掛節(jié)點的問題.

當(dāng)四叉樹網(wǎng)格嚴(yán)格遵循2:1 的規(guī)則時,即相鄰網(wǎng)格的尺寸比最大不能超過2:1,稱該種網(wǎng)格為平衡四叉樹網(wǎng)格.圖5 為一簡單平衡四叉樹網(wǎng)格劃分示意圖,劃分網(wǎng)格的基本思想是將二維結(jié)構(gòu)離散為若干個正方形子網(wǎng)格,再對子網(wǎng)格重復(fù)遞歸離散,達(dá)到尺寸要求后方可終止離散.該網(wǎng)格離散算法能實現(xiàn)不同網(wǎng)格尺寸層級的快速劃分及平滑過渡,從米級到毫米級的網(wǎng)格僅需10 次遞歸離散(210=1024).對于各向同性材料,采用平衡四叉樹網(wǎng)格時,只存在圖6 中的6 種單元類型.根據(jù)前述,在損傷計算過程中,SBFEM 子域的特征值、位移模態(tài)、應(yīng)力模態(tài)是不變的,且每一荷載步子域的剛度矩陣都可從初始剛度矩陣計算得到,只有積分常數(shù)c需要根據(jù)每一步的位移場進(jìn)行計算更新.將上述6 種單元類型的單元計算信息(特征值、位移模態(tài)、初始剛度矩陣等)預(yù)先計算并存儲,在計算時可直接調(diào)用,有利于節(jié)省計算時間.故將SBFEM 與四叉樹網(wǎng)格結(jié)合進(jìn)行數(shù)值計算時,不僅使前處理階段變得方便省時,還有效地提高了計算效率.

圖5 平衡四叉樹網(wǎng)格劃分示意圖Fig.5 Balanced quadtree discretization process

圖6 平衡四叉樹網(wǎng)格的6 種單元類型Fig.6 The six cell types of the balanced quadtree mesh

4 數(shù)值算例

本文通過兩個典型試件的準(zhǔn)靜態(tài)開裂模擬,并與試驗結(jié)果和其他文獻(xiàn)結(jié)果進(jìn)行比較,以此驗證本文方法的正確性與有效性.為了充分說明本文方法在計算效率上的提升,使用配備Intel (R) Core (TM)i5-4200 H CPU @ 2.80 GHz 的計算機(jī)進(jìn)行計算,該處理器的性能與文獻(xiàn)[15] 所采用的Intel (R) Core(TM) i5-4210 M CPU @ 2.60 GH 接近.對于大部分模型,損傷只會出現(xiàn)在易損的局部區(qū)域內(nèi),同時也為了保證式(17)的數(shù)值積分精度,因此只需在損傷過程區(qū)劃分足夠精細(xì)的網(wǎng)格,其余部分的網(wǎng)格滿足精度要求即可.所有算例均考慮為平面應(yīng)力問題.

4.1 預(yù)制缺口三點彎曲梁

圖7 所示為一預(yù)制缺口三點彎曲梁,其幾何尺寸、邊界條件和荷載信息如圖所示.該算例屬于張開型失效破壞,Rots[32]對該梁進(jìn)行了相關(guān)試驗研究.材料參數(shù)根據(jù)文獻(xiàn)[32]選取,彈性模量E=20 GPa,泊松比v=0.2.非局部宏微觀損傷模型中的參數(shù)根據(jù)文獻(xiàn)[15]選取,影響域半徑l=5 mm,臨界伸長率λc=1.2 × 10-4,γ=1000,p=3,q=4.劃分網(wǎng)格時,對缺口附近一定范圍內(nèi)的區(qū)域進(jìn)行網(wǎng)格加密處理.為了研究網(wǎng)格敏感性問題,采用如圖8 所示3 種不同疏密程度的四叉樹網(wǎng)格進(jìn)行計算.3 種網(wǎng)格在損傷過程區(qū)的網(wǎng)格尺寸均小于l/5,表1 給出了3 種網(wǎng)格的具體信息.

表1 預(yù)制缺口三點彎曲梁網(wǎng)格信息Table 1 Mesh information of the notched beam

圖7 預(yù)制缺口三點彎曲梁(單位:mm)Fig.7 Three-point bending of the notched beam (unit:mm)

圖8 預(yù)制缺口三點彎曲梁四叉樹網(wǎng)格Fig.8 Quadtree mesh of a notched beam

3 種不同網(wǎng)格計算得到的最終損傷分布如圖9(a)～圖9(c)所示,可以看出不同網(wǎng)格的裂紋形態(tài)與走向基本相同.與文獻(xiàn)[25]中采用的積分型非局部損傷模型所計算的損傷分布(圖9(d))相比,本文的損傷分布是一條細(xì)窄帶,能更形象表征裂紋形態(tài),而文獻(xiàn)中的損傷帶明顯較寬,對于裂紋形態(tài)的表征不夠準(zhǔn)確.與文獻(xiàn)[33]中采用內(nèi)聚力模型的計算結(jié)果(圖9(e))相比,文獻(xiàn)中的裂紋走向會受到網(wǎng)格影響,其主要原因是內(nèi)聚單元是插入到單元邊界的,裂紋擴(kuò)展只能沿單元邊界擴(kuò)展,故內(nèi)聚力模型的裂紋路徑易受網(wǎng)格影響,表現(xiàn)出網(wǎng)格敏感性,而本文方法的裂紋路徑基本不受網(wǎng)格影響.圖10 給出預(yù)制缺口三點彎曲梁的荷載位移曲線,可以看出三種網(wǎng)格的荷載位移曲線具有良好的一致性.可見在網(wǎng)格尺寸小于l/5 時,本文方法能有效克服經(jīng)典局部損傷力學(xué)中的網(wǎng)格敏感性問題.圖10 還給出了與試驗結(jié)果、文獻(xiàn)計算結(jié)果的比較,本文荷載位移曲線與文獻(xiàn)[32]的試驗曲線范圍基本吻合.在與其他文獻(xiàn)的對比中,文獻(xiàn)[25]采用的是基于SBFEM 的積分型非局部損傷模型,文獻(xiàn)[15]采用的是基于有限元的非局部宏微觀損傷模型,文獻(xiàn)[18]采用的是統(tǒng)一相場模型.可以看出本文方法與其他文獻(xiàn)中的方法都能較準(zhǔn)確地捕捉荷載位移曲線,只是在下降段存在些許差別.

圖9 預(yù)制缺口三點彎曲梁最終損傷分布云圖Fig.9 Final damage distribution of three-point bending failure of the notched beam

圖10 預(yù)制缺口三點彎曲梁荷載位移曲線Fig.10 Load-displacement curve for different mesh sizes

圖11 給出了4 個典型階段的損傷與位移分布云圖,分別對應(yīng)于圖10 網(wǎng)格A 曲線上所標(biāo)記出的四個階段.從四個階段的損傷分布云圖可以看出,裂紋從缺口尖端開始垂直向梁頂部擴(kuò)展,與在試驗中觀察到的現(xiàn)象一致.當(dāng)達(dá)到峰值荷載(階段Ⅱ)時,已經(jīng)能看出較為明顯的損傷,即開始出現(xiàn)肉眼可見的裂紋.在下降段時期(階段Ⅱ之后),裂紋進(jìn)入不穩(wěn)定的發(fā)展階段,裂紋兩側(cè)的位移場也表現(xiàn)出明顯的不連續(xù)性.

圖11 預(yù)制缺口三點彎曲梁典型階段損傷與位移分布云圖 (mesh A)Fig.11 Damage and displacement field of the notched beam (mesh A)

圖12 給出了網(wǎng)格A 裂紋附近的3 個物質(zhì)點M(219.69,55.31),N(222.19,55.31),O(224.69,55.31)所連物質(zhì)鍵在4 個典型階段的損傷情況.從圖中可以看出垂直于裂紋方向的物質(zhì)鍵最先出現(xiàn)損傷,且損傷程度更大.3 個物質(zhì)點是處于同一水平線上的,距離裂紋越近的點,受損的物質(zhì)鍵越多.圖13是3 個物質(zhì)點的拓?fù)鋼p傷發(fā)展曲線,可以看出M和N兩點的拓?fù)鋼p傷在階段Ⅲ之后已趨于穩(wěn)定,且損傷程度較小.是因為M和N兩點只有少量垂直于裂紋方向的物質(zhì)鍵受損,且階段Ⅲ后這些物質(zhì)鍵的損傷程度已達(dá)極限.而O點的拓?fù)鋼p傷在階段Ⅳ后才趨于穩(wěn)定,且損傷程度已趨于完全損傷.其主要原因是O點靠近裂紋,附近的位移場具有明顯的不連續(xù)性,有更多的物質(zhì)鍵受到更大程度的損傷.從圖中還可以看出,距裂紋較遠(yuǎn)的點(如M點)的拓?fù)鋼p傷結(jié)果僅取決于部分方向上的物質(zhì)鍵的損傷程度(大部分方向的物質(zhì)鍵的損傷程度很小),計算的損傷值相對較小,從完全損傷到較小損傷過渡非常快.而常規(guī)非局部損傷模型損傷計算通過非局部變量(如等效應(yīng)變或塑性應(yīng)變)進(jìn)行,這些非局部變量計算與所有方向上的點有關(guān),導(dǎo)致非局部變量的值偏大,最終導(dǎo)致計算的損傷值偏大,計算的損傷帶就較寬.而本文采用近場動力學(xué)的物質(zhì)鍵失效的思路的損傷計算有別于常規(guī)有限元的損傷計算.同時,本文的損傷模型引進(jìn)了統(tǒng)一相場的能量退化函數(shù),使得從完好到完全損傷的過渡更為迅速.最終體現(xiàn)為本文方法的損傷帶是一條細(xì)窄帶.

圖12 點M(左)、N(中)、O(右)在典型階段的鍵損傷分布Fig.12 Damage of bonds of points M (left),N(middle) and O(right) in typical stage

圖12 點M(左)、N(中)、O(右)在典型階段的鍵損傷分布(續(xù))Fig.12 Damage of bonds of points M (left),N(middle) and O(right) in typical stage (continued)

圖13 點M,N,O 的拓?fù)鋼p傷演化Fig.13 Topological damage evolution of points M,N and O

表2 給出了本文與文獻(xiàn)[15]的單元數(shù)和計算時間對比,其中文獻(xiàn)[15]采用有限元與三角形網(wǎng)格進(jìn)行計算.因本文模型的網(wǎng)格A、網(wǎng)格B 和網(wǎng)格C 的計算結(jié)果幾乎完全相同,以網(wǎng)格A 與文獻(xiàn)[15]的網(wǎng)格A 比較如下,本文的單元數(shù)為3530 個,文獻(xiàn)[15]的單元數(shù)為5712 個,這些單元主要集中在損傷過程區(qū),前者只有后者的62%.從表中可以看出,在處理器性能與子域(單元)數(shù)相當(dāng)?shù)臈l件下,本文的計算時間較文獻(xiàn)[15]大幅縮短.說明相較于常規(guī)有限元而言,將比例邊界有限元與四叉樹網(wǎng)格結(jié)合應(yīng)用于非局部宏微觀損傷模擬,計算效率可得到有效提高.其主要原因是四叉樹網(wǎng)格可直接與比例邊界有限元結(jié)合使用,而無需處理懸掛節(jié)點問題,且六種單元的計算信息可預(yù)先計算并存儲,計算時直接調(diào)用即可,極大地提高了計算效率.

表2 本文與文獻(xiàn)[15]的單元數(shù)和計算時間Table 2 Number of elements and calculation time in this paper are compared with Ref.[15]

4.2 L 形混凝土板

圖14 為一個L 形混凝土板,其幾何尺寸、邊界條件和荷載信息如圖所示.該算例屬于混合型失效破壞,且沒有預(yù)設(shè)初始缺陷,Winkler等[34]對該L 形混凝土板進(jìn)行了相關(guān)試驗研究.材料參數(shù)根據(jù)文獻(xiàn)[34]選取,彈性模量E=25.85 GPa,泊松比v=0.18.非局部宏微觀損傷模型中的參數(shù)根據(jù)文獻(xiàn)[15]選取,影響域半徑l=10 mm,臨界伸長率λc=1.05 × 10-4,γ=850,p=4,q=4.劃分網(wǎng)格時,對折角處至左自由邊界一定范圍內(nèi)的區(qū)域進(jìn)行網(wǎng)格加密處理.為了研究網(wǎng)格敏感性問題,采用如圖15 所示3 種不同疏密程度的四叉樹網(wǎng)格進(jìn)行計算.3 種網(wǎng)格在損傷過程區(qū)的網(wǎng)格尺寸均小于l/5,表3 給出了3 種網(wǎng)格的具體信息.

圖14 L 形混凝土板(單位:mm)Fig.14 Winkler L-shaped plane (unit:mm)

表3 L 形混凝土板網(wǎng)格信息Table 3 Mesh information of L-shaped plane

圖15 L 形混凝土板四叉樹網(wǎng)格Fig.15 Quadtree mesh of L-shaped plane

圖16 給出了L 形混凝土板最終損傷分布與試驗裂紋擴(kuò)展路徑,可以看出不同網(wǎng)格所計算的裂紋形態(tài)與走向是基本相同的,并且與試驗裂紋路徑吻合.圖17 給出了L 形混凝土板的荷載位移曲線,可以看出三種網(wǎng)格的荷載位移曲線基本一致.再次說明當(dāng)網(wǎng)格尺寸小于l/5 時,本文方法不存在網(wǎng)格敏感性問題.從圖17 還可看出,本文荷載位移曲線基本處于文獻(xiàn)[34]的試驗曲線范圍內(nèi),Winkler等[34]同時還采用基于有限元的各向同性損傷模型對裂縫開裂進(jìn)行了數(shù)值模擬(圖中黑色虛線).文獻(xiàn)[15]采用的是基于有限元的非局部宏微觀損傷模型,文獻(xiàn)[18]采用的是統(tǒng)一相場模型.可以看出本文曲線與文獻(xiàn)[15,34] 所計算得到的曲線差距很小,文獻(xiàn)[18]的曲線雖然與試驗曲線差距偏大,但總體趨勢是正確的.

圖16 L 形混凝土板最終損傷分布與試驗裂紋路徑Fig.16 Crack patterns (damage) for different mesh sizes

圖16 L 形混凝土板最終損傷分布與試驗裂紋路徑(續(xù))Fig.16 Crack patterns (damage) for different mesh sizes (continued)

圖17 L 形混凝土板荷載-位移曲線Fig.17 Load-displacement curve for different mesh sizes

圖18 給出了四個典型階段的損傷分布與位移分布云圖,分別對應(yīng)于圖17 網(wǎng)格A 曲線上所標(biāo)記出的四個階段.從四個階段的損傷分布云圖可以看出,裂紋起始于折角點,在向左上方擴(kuò)展一段距離后,便水平向左自由邊界擴(kuò)展,與在試驗中觀察到的現(xiàn)象一致.類似于前述三點彎曲梁觀察到的現(xiàn)象,L 形板在達(dá)到峰值荷載(階段Ⅱ)時,從圖18(c)～圖18(d)已能觀察到較為明顯的裂紋.在下降段時期(階段Ⅱ之后),裂紋兩側(cè)的位移場也表現(xiàn)出明顯的不連續(xù)性,說明裂紋已進(jìn)入不穩(wěn)定的發(fā)展階段.

圖18 L 形混凝土板典型階段損傷與位移分布云圖 (mesh A)Fig.18 Damage and displacement field of the L-shaped plane (mesh A)

圖18 L 形混凝土板典型階段損傷與位移分布云圖 (mesh A)(續(xù))Fig.18 Damage and displacement field of the L-shaped plane (mesh A)(continued)

表4 給出了本文與文獻(xiàn)[15]的計算單元數(shù)和計算時間的比較,其中文獻(xiàn)[15]采用三角形單元進(jìn)行計算.與前類似,選擇本文模型的網(wǎng)格A 與文獻(xiàn)[15]的網(wǎng)格A 進(jìn)行比較如下,前者單元數(shù)為6396 個,后者單元數(shù)為10 902 個,前者大約為后者的59%.可以明顯看出,本文的計算時間明顯少于文獻(xiàn)[15],只有后者的19%.再次說明將比例邊界有限元與四叉樹網(wǎng)格結(jié)合使用,對于提升非局部宏微觀損傷模擬的計算效率是非常顯著的.

表4 本文與文獻(xiàn)[15]的單元數(shù)和計算時間Table 4 Number of element and calculation time in this paper are compared with Ref.[15]

圖19 給出了網(wǎng)格A 裂紋附近的3 個物質(zhì)點P(235.35,247.07),Q(238.35,250.98),R(235.35,254.88)所連物質(zhì)鍵在4 個典型階段的損傷情況.圖20給出了這3 個物質(zhì)點的拓?fù)鋼p傷發(fā)展曲線.與前述算例類似,垂直于裂紋方向的物質(zhì)鍵最先出現(xiàn)損傷,且越靠近裂紋的物質(zhì)點所連物質(zhì)鍵受損程度越大,受損的物質(zhì)鍵數(shù)量越多,則相應(yīng)的拓?fù)鋼p傷也越大.且大程度的拓?fù)鋼p傷集中于位移場不連續(xù)的位置,因而損傷帶表現(xiàn)為一條細(xì)窄帶.

圖19 點P(左)、Q(中)、R(右)在典型階段的鍵損傷分布Fig.19 Damage of bonds of points P (left),Q(middle) and R(right) in typical stage

圖19 點P(左)、Q(中)、R(右)在典型階段的鍵損傷分布(續(xù))Fig.19 Damage of bonds of points P (left),Q(middle) and R(right) in typical stage (continued)

圖20 點P,Q,R 的拓?fù)鋼p傷演化Fig.20 Topological damage evolution of points P,Q and R

5 結(jié)論

本文將比例邊界有限元與非局部宏微觀損傷模型結(jié)合,采用四叉樹離散技術(shù)進(jìn)行網(wǎng)格劃分,并將其應(yīng)用于準(zhǔn)脆性材料的裂紋模擬中.通過兩個典型算例的驗證分析,得到如下主要結(jié)論.

(1) 本文方法適用于各向同性準(zhǔn)脆性材料的準(zhǔn)靜態(tài)開裂模擬,能準(zhǔn)確計算裂紋擴(kuò)展路徑與荷載-變形曲線.對于沒有預(yù)制初始缺陷的問題,裂紋的起裂與擴(kuò)展也能自發(fā)進(jìn)行.

(2) 本文模型是一種引入非局部思想的非局部化模型,當(dāng)損傷過程區(qū)的網(wǎng)格尺寸小于影響域半徑的1/5 時,不存在經(jīng)典局部損傷力學(xué)中的網(wǎng)格敏感性問題.

(3) 比例邊界有限元與四叉樹網(wǎng)格結(jié)合,不存在懸掛節(jié)點問題,實現(xiàn)了快速、高質(zhì)量的多級網(wǎng)格劃分與尺寸過渡.六種單元類型的計算信息可預(yù)先計算并存儲,有效地提高了計算效率.

本文的研究工作已初步驗證了比例邊界有限元與四叉樹網(wǎng)格在非局部宏微觀損傷模擬中的可行性與高效性.但還存在一系列值得深入研究的問題需要在未來解決,例如非均質(zhì)材料的損傷模擬、各向異性損傷模擬、動態(tài)損傷模擬等問題.