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

杜成斌 黃文倉 江守燕

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

引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

則有

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

將Hamilton 矩陣進行特征值分解,對于每個子域都滿足下式關系

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

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

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

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

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

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

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

2 非局部宏微觀損傷模型及其在SBFEM中的實現

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

2.1 微細觀損傷與拓撲損傷

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

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

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

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

式中,γ 為材料參數,與材料的微觀結構有關,可通過標準試件試驗確定.并且 γ 表征損傷發展的速度,γ的取值越大,損傷發展越快,脆性性質越明顯,圖3給出了不同取值的 γ 對應的微細觀損傷演化曲線.加載歷史參數 κ 為物質鍵的歷史最大超越伸長率,其定義為[15-16]

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

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

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

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

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

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

2.2 能量退化函數

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

式中,p和q均為退化參數,其取值刻畫了材料的宏觀脆性程度,可以認為是材料的物性參數,可以通過標準試件試驗確定[16,30].

2.3 非局部宏微觀損傷模型在SBFEM 中的實現

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

其中,σ 為應力張量,ε 為應變張量.在非局部宏微觀損傷模型中,通過能量退化函數g=g(Ω) 的引入,建立起拓撲損傷與應力應變之間的聯系,則本構關系轉變為

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

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

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

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

3 四叉樹網格

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

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

圖5 平衡四叉樹網格劃分示意圖Fig.5 Balanced quadtree discretization process

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

4 數值算例

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

4.1 預制缺口三點彎曲梁

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

表1 預制缺口三點彎曲梁網格信息Table 1 Mesh information of the notched beam

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

圖8 預制缺口三點彎曲梁四叉樹網格Fig.8 Quadtree mesh of a notched beam

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

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

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

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

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

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

圖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(右)在典型階段的鍵損傷分布(續)Fig.12 Damage of bonds of points M (left),N(middle) and O(right) in typical stage (continued)

圖13 點M,N,O 的拓撲損傷演化Fig.13 Topological damage evolution of points M,N and O

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

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

4.2 L 形混凝土板

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

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

表3 L 形混凝土板網格信息Table 3 Mesh information of L-shaped plane

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

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

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

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

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

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

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

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

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

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

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

圖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(右)在典型階段的鍵損傷分布(續)Fig.19 Damage of bonds of points P (left),Q(middle) and R(right) in typical stage (continued)

圖20 點P,Q,R 的拓撲損傷演化Fig.20 Topological damage evolution of points P,Q and R

5 結論

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

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

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

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

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