基于碎點(diǎn)法的動(dòng)態(tài)斷裂分析1)

沈?qū)毈?王松 李明凈 董雷霆

(北京航空航天大學(xué)航空科學(xué)與工程學(xué)院,北京 100191)

引言

沖擊防護(hù)結(jié)構(gòu)廣泛存在于航空航天、汽車、船舶、核能乃至國防等工程領(lǐng)域.這類結(jié)構(gòu)在抵抗如撞擊、爆炸等極端的沖擊載荷時(shí)可能會發(fā)生動(dòng)態(tài)斷裂直至最終破壞.準(zhǔn)確預(yù)測結(jié)構(gòu)在動(dòng)態(tài)載荷下的斷裂破壞行為對防護(hù)結(jié)構(gòu)的設(shè)計(jì)和改進(jìn)具有重要意義.數(shù)值仿真方法是預(yù)測結(jié)構(gòu)動(dòng)態(tài)斷裂的重要手段,學(xué)者在這個(gè)方面已經(jīng)進(jìn)行了長期、大量的研究,然而由于動(dòng)態(tài)斷裂問題的復(fù)雜性,這仍然是計(jì)算固體力學(xué)領(lǐng)域一個(gè)具有挑戰(zhàn)性的研究方向.

有限元法是當(dāng)前工程領(lǐng)域應(yīng)用最廣泛的數(shù)值仿真方法,有限元法基于伽遼金弱形式且使用可以使用簡單多項(xiàng)式的形函數(shù),具有穩(wěn)定性強(qiáng)、數(shù)值積分簡單等優(yōu)點(diǎn)[1].然而有限元法在模擬動(dòng)態(tài)斷裂導(dǎo)致的結(jié)構(gòu)破壞時(shí)存在兩個(gè)方面的難題.第一個(gè)方面,有限元法本質(zhì)上是基于位移連續(xù)性假設(shè)的,因此難以在模型中顯式引入裂紋.當(dāng)前解決這個(gè)問題的一種方法是使用網(wǎng)格刪除技術(shù),然而刪除網(wǎng)格會導(dǎo)致和物理過程不一致的材料損失[2].另一種方法是使用內(nèi)聚力單元模擬內(nèi)界面失效,但是固有型內(nèi)聚力單元存在人工柔度問題,非固有型內(nèi)聚力單元需要使用較大的內(nèi)界面懲罰系數(shù),在模擬復(fù)雜斷裂問題時(shí)存在困難[3-4].擴(kuò)展有限元法能在有限元模型中顯式引入裂紋,是一種受到廣泛關(guān)注的計(jì)算方法,但該方法在引入裂紋時(shí)會引入額外的自由度,且通常需要對模型中的裂紋區(qū)域進(jìn)行網(wǎng)格細(xì)化和重構(gòu),計(jì)算效率和穩(wěn)定性有待進(jìn)一步提高[5].第二個(gè)方面,使用有限元法模擬復(fù)雜動(dòng)態(tài)斷裂問題時(shí)往往會伴隨裂紋附近局部區(qū)域的網(wǎng)格畸變,由此導(dǎo)致顯式動(dòng)力學(xué)求解的臨界時(shí)間步長顯著減小,解的精度下降[6-9].過高的網(wǎng)格畸變甚至?xí)蜇?fù)體積等問題導(dǎo)致計(jì)算終止[9].當(dāng)前通常的解決辦法是在計(jì)算過程中刪除模型中的畸變,但這同樣會導(dǎo)致與物理過程不符的材料損失.

以上所述的傳統(tǒng)有限元方法在模擬動(dòng)態(tài)斷裂中面臨的諸多制約,其根源在于有限元方法使用單元離散求解域,并基于單元節(jié)點(diǎn)計(jì)算形函數(shù),因此有限元對網(wǎng)格具有很強(qiáng)的依賴性.無網(wǎng)格方法是部分或完全擺脫網(wǎng)格依賴的一類數(shù)值方法,其核心思想是采用節(jié)點(diǎn)群對結(jié)構(gòu)進(jìn)行離散,并使用節(jié)點(diǎn)描述被離散結(jié)構(gòu)的力學(xué)響應(yīng)[10].因?yàn)闊o網(wǎng)格類方法能有效避免傳統(tǒng)有限元法普遍存在的網(wǎng)格畸變問題,因此這類方法一直是國內(nèi)外計(jì)算力學(xué)界的研究熱點(diǎn)之一.從20世紀(jì)90年代起,學(xué)者們提出的無網(wǎng)格方法已達(dá)20余種,這些無網(wǎng)格法大體上可以分為兩種:粒子類方法和基于弱形式的無網(wǎng)格方法.

光滑粒子流體動(dòng)力學(xué)(SPH)方法[11-12]是沖擊領(lǐng)域應(yīng)用較為廣泛的一種粒子類無網(wǎng)格法,該方法利用光滑核函數(shù)近似模擬場函數(shù),其離散控制方程是配點(diǎn)形式的[13].每個(gè)SPH 粒子具有質(zhì)量、速度、體積等物理屬性,這種特性使得SPH 在顯式動(dòng)力學(xué)仿真中具有很大優(yōu)勢.但由于SPH 基于配點(diǎn)的強(qiáng)形式,其穩(wěn)定性難以獲得嚴(yán)格的數(shù)學(xué)證明,因此SPH 往往需要較大的支持域,導(dǎo)致SPH 法的計(jì)算量通常高于有限元法[14].物質(zhì)點(diǎn)法(MPM)是另一種代表性的粒子類方法.MPM 采用歐拉和拉格朗日雙重描述將物質(zhì)離散為一組在空間網(wǎng)格上運(yùn)動(dòng)的質(zhì)點(diǎn),這些質(zhì)點(diǎn)同樣攜帶質(zhì)量、體積、速度等物理信息.動(dòng)量方程在歐拉描述的空間網(wǎng)格中離散,避免了大變形造成拉格朗日網(wǎng)格畸變的問題[15].文獻(xiàn)[16-18]采用物質(zhì)點(diǎn)法在物體超高速碰撞、沖擊侵徹、爆炸、動(dòng)態(tài)斷裂等方面開展了一系列研究,取得了重要進(jìn)展.

Atluti等[19-20]提出的無網(wǎng)格局部彼得羅夫-伽遼金(MLPG)方法和Belytschko等[21-22]提出的無單元伽遼金(EFG)方法是兩種較為典型的弱形式無網(wǎng)格方法.這兩種無網(wǎng)格法基于弱形式方程,因此具有較好的數(shù)值穩(wěn)定性[19,21,23].然而該類方法的形函數(shù)通常為復(fù)雜有理函數(shù),這導(dǎo)致了在積分時(shí)無法籍由簡單的高斯積分直接計(jì)算,積分方法不當(dāng)甚至?xí)绊懛椒ǖ姆€(wěn)定性.文獻(xiàn)[24-26]提出一系列積分方法,在降低計(jì)算量、提高積分精度及方法穩(wěn)定性等方面具有明顯效果.

近年來,針對動(dòng)態(tài)載荷作用下的復(fù)雜結(jié)構(gòu)破壞問題,國內(nèi)外學(xué)者在近場動(dòng)力學(xué)、邊界元、間斷伽遼金有限元法等方面也開展了一定工作.劉立勝等采用近場動(dòng)力學(xué)開展了濕熱環(huán)境下復(fù)合材料沖擊損傷及動(dòng)態(tài)斷裂的研究[27],章青等進(jìn)行了爆炸載荷作用下混凝土結(jié)構(gòu)破壞過程的近場動(dòng)力學(xué)模擬[28-29],Han等[30]研究了耦合損傷力學(xué)與近場動(dòng)力學(xué)的材料失效仿真方法,高效偉等[31]采用邊界元法分析功能梯度材料動(dòng)態(tài)斷裂力學(xué)問題,于福臨等開展了基于間斷伽遼金方法的船體板架爆炸沖擊響應(yīng)數(shù)值模擬研究[32],相關(guān)的研究工作仍在進(jìn)行當(dāng)中.

綜上所述,國內(nèi)外學(xué)者針對動(dòng)態(tài)斷裂問題已經(jīng)開展了多種數(shù)值仿真方法的研究,但由于該問題物理過程的復(fù)雜性,對動(dòng)態(tài)斷裂過程的準(zhǔn)確仿真預(yù)測仍然是一項(xiàng)具有挑戰(zhàn)性的工作.

近年來,作者所在研究團(tuán)隊(duì)提出了一種不連續(xù)型無網(wǎng)格方法:碎點(diǎn)法(FPM)[33],該方法一方面參考弱形式無網(wǎng)格方法,采用弱形式方程,具有較好的數(shù)值穩(wěn)定性,且基于支持域節(jié)點(diǎn)群定義形函數(shù),使子域具有抵抗畸變的能力;另一方面參考間斷伽遼金有限元法,不對相鄰子域間的位移作連續(xù)性要求,因此該方法易于在模型中顯式引入裂紋.

碎點(diǎn)法使用空間中的節(jié)點(diǎn)群對問題域進(jìn)行離散,基于節(jié)點(diǎn)群劃分具有明確幾何形狀的子域,如圖1(a)所示.碎點(diǎn)法模型中僅考慮節(jié)點(diǎn)的自由度,子域頂點(diǎn)信息由節(jié)點(diǎn)插值獲得.節(jié)點(diǎn)反映的是相應(yīng)子域的物理狀態(tài),因此節(jié)點(diǎn)具有體積、質(zhì)量、速度等明確的物理含義.對于碎點(diǎn)法模型中的任意一個(gè)節(jié)點(diǎn)P0及其子域E0,其支持域?yàn)樽佑駿0的所有相鄰子域E1～E6,則子域E0中的位移形函數(shù)是基于支持域節(jié)點(diǎn)群P1～P6定義的,如圖1(b)所示.由于碎點(diǎn)法使用支持域節(jié)點(diǎn)群而非子域頂點(diǎn)構(gòu)建形函數(shù),計(jì)算過程中子域形狀的畸變不會對形函數(shù)的計(jì)算造成影響,因此碎點(diǎn)法具有抵抗子域形狀畸變的能力.當(dāng)然碎點(diǎn)法的支持域定義具有很大的靈活性,可以根據(jù)實(shí)際計(jì)算的需要進(jìn)行修改,本文僅以圖1(b)所示的方法為例進(jìn)行說明.

圖1 (a) 碎點(diǎn)法模型的離散和(b) 支持域及其節(jié)點(diǎn)群的定義Fig.1 (a) Discretization of FPM model and(b) definition of support range and its point cloud

碎點(diǎn)法不對相鄰子域間的位移形函數(shù)和試函數(shù)作連續(xù)性要求,即所使用的形函數(shù)和試函數(shù)可以是分片連續(xù)的,如圖2 所示.由于只要求函數(shù)在子域中局部連續(xù),碎點(diǎn)法可以使用簡單多項(xiàng)式形式的位移形函數(shù)和試函數(shù),因此子域內(nèi)的積分計(jì)算可以使用簡單的高斯積分進(jìn)行求解,從而避免了大多數(shù)弱形式無網(wǎng)格法因需要使用復(fù)雜有理數(shù)形函數(shù)而導(dǎo)致積分困難的問題.另一方面,由于碎點(diǎn)法使用的位移試函數(shù)是分片連續(xù)的,可以自然地在任意相鄰子域之間的內(nèi)部界面處顯式引入裂紋,因此碎點(diǎn)法具備模擬大范圍復(fù)雜裂紋萌生和擴(kuò)展問題的能力.在碎點(diǎn)法模型中顯式引入裂紋的具體方法將在下文中詳細(xì)介紹.

圖2 (a) 碎點(diǎn)法形函數(shù)示意圖(b) 位移場為的碎點(diǎn)法位移試函數(shù)[34]Fig.2(a) Shapefunctionof FPM(b) FPM’strialfunction of a field of displacement:

碎點(diǎn)法的方程是基于伽遼金弱形式構(gòu)建的,因此該方法具有較好的數(shù)值穩(wěn)定性.然而,碎點(diǎn)法的位移形函數(shù)和試函數(shù)是分片連續(xù)的,模型中相鄰子域間的位移連續(xù)性缺乏約束,因此方法的一致性無法得到保證.為了解決這個(gè)問題,參考間斷伽遼金有限元法,在弱形式方程中任意內(nèi)界面處引入數(shù)值通量修正,以弱形式約束相鄰子域的位移連續(xù)性.在作者所在團(tuán)隊(duì)的前期工作中,已經(jīng)證明了在弱形式方程中引入內(nèi)界面數(shù)值通量修正可以有效保證碎點(diǎn)法計(jì)算的一致性[33].

綜上所述,碎點(diǎn)法具有如下優(yōu)勢:(1)采用基于節(jié)位移形函數(shù),因此具有抵抗畸變的能力;(2)放棄了相鄰子域間的形函數(shù)連續(xù)性要求,因此可以采用簡單多項(xiàng)式形式的位移試函數(shù)并用高斯積分求解弱形式方程,且易于在模型中顯式引入裂紋;(3)子域具有質(zhì)量、體積、速度等物理屬性,因此便于顯式動(dòng)力學(xué)的計(jì)算.基于上述特點(diǎn),作者認(rèn)為碎點(diǎn)法適合用于預(yù)測結(jié)構(gòu)的動(dòng)態(tài)斷裂行為.

碎點(diǎn)法已經(jīng)在多個(gè)領(lǐng)域得到了檢驗(yàn)和應(yīng)用.Yang等[34]應(yīng)用碎點(diǎn)法進(jìn)行準(zhǔn)靜態(tài)應(yīng)力分析和斷裂分析,驗(yàn)證了碎點(diǎn)法求解上述問題的有效性和準(zhǔn)確性;Wang等[35]應(yīng)用碎點(diǎn)法預(yù)測含U 型缺口脆性試件的斷裂強(qiáng)度和裂紋形貌,證明碎點(diǎn)法可以有效分析復(fù)雜的準(zhǔn)靜態(tài)脆性斷裂問題;Guan等[36-37]證明了在分析撓曲電材料的斷裂時(shí),碎點(diǎn)法相比于傳統(tǒng)有限元法具有計(jì)算簡單、易于顯式引入裂紋的特點(diǎn);Guan等[38-39]應(yīng)用碎點(diǎn)法求解非均質(zhì)材料與復(fù)雜薄壁結(jié)構(gòu)的瞬態(tài)熱傳導(dǎo)問題,展示了碎點(diǎn)法在解決復(fù)雜問題域及邊界的瞬態(tài)熱傳導(dǎo)問題方面的準(zhǔn)確性、高效性及穩(wěn)定性.

本文旨在提出基于碎點(diǎn)法核心思想的動(dòng)力學(xué)碎點(diǎn)法理論并開發(fā)相應(yīng)程序,同時(shí)驗(yàn)證該方法分析應(yīng)力波傳播和動(dòng)態(tài)裂紋擴(kuò)展的有效性.本文的理論推導(dǎo)和算例驗(yàn)證都是基于二維空間小變形假設(shè).本文的結(jié)構(gòu)安排如下:第一章介紹碎點(diǎn)法的基本理論,包括碎點(diǎn)法基本理論、求解域離散方法、弱形式動(dòng)量方程及其離散、顯式動(dòng)力學(xué)求解格式;第二章通過算例驗(yàn)證動(dòng)力學(xué)碎點(diǎn)法在預(yù)測應(yīng)力波傳播和動(dòng)態(tài)裂紋擴(kuò)展方面的能力;第三章對本文做出結(jié)論.

1 碎點(diǎn)法基本理論

1.1 求解域離散

碎點(diǎn)法采用布置節(jié)點(diǎn)的方式對求解域 Ω 進(jìn)行離散,在建立碎點(diǎn)法方程時(shí)僅考慮節(jié)點(diǎn)處的自由度.基于節(jié)點(diǎn)進(jìn)一步將求解域 Ω 劃分若干個(gè)子域 Ωi,每個(gè)子域內(nèi)部包含一個(gè)節(jié)點(diǎn),子域具有明確的幾何形狀,且相鄰子域之間互不相交重疊.碎點(diǎn)法的離散方式如圖1(a)所示.

碎點(diǎn)法的離散有兩種基本方式,其一是先在求解域中布置節(jié)點(diǎn),然后基于節(jié)點(diǎn)定義Voronoi 多邊形作為子域,如圖1(a)所示.另一種是借助成熟的有限元前處理軟件,先在求解域中劃分有限元網(wǎng)格,以有限元單元作為碎點(diǎn)法的子域,然后在子域中定義節(jié)點(diǎn).

對于碎點(diǎn)法模型中的任意子域,在子域節(jié)點(diǎn)處對其位移函數(shù)u進(jìn)行泰勒級數(shù)展開,從而將該子域內(nèi)的位移函數(shù)表示為節(jié)點(diǎn)位移和節(jié)點(diǎn)位移梯度的函數(shù),由此便定義了子域內(nèi)的多項(xiàng)式形式的位移試函數(shù),如式(1)所示.以包含節(jié)點(diǎn)P0的子域E0為例,在P0處進(jìn)行泰勒級數(shù)展開,則子域E0內(nèi)的位移試函數(shù)如下

節(jié)點(diǎn)處的位移梯度可以采用多種方法進(jìn)行求解,本文所采用的是廣義有限差分法.首先定義目標(biāo)節(jié)點(diǎn)P0和子域E0的支持域:本文定義與子域E0相鄰且具有公共邊界的子域作為支持域,這些子域的節(jié)點(diǎn)P1,P2,···,Pm組成支持域節(jié)點(diǎn)群,如圖1(b) 所示.在定義了子域的支持域之后,可通過廣義有限差分法來求解位移梯度.

定義L2范數(shù)J如下式所示

其中a為位移梯度向量,u0和um分別為組裝后的P0和P1～Pm節(jié)點(diǎn)位移向量,A是包含P0和P1～Pm之間坐標(biāo)差的張量,W是權(quán)函數(shù)張量,具體表達(dá)式如下

其中,uE為包含P0節(jié)點(diǎn)支持域節(jié)點(diǎn)群(P0～Pm)位移的組裝位移向量,I1和I2為兩個(gè)轉(zhuǎn)換張量,表達(dá)式如下

將式(4)代入式(3)可得節(jié)點(diǎn)P0處位移梯度向量a與支持域節(jié)點(diǎn)群位移向量uE的表達(dá)式如下

將式(5)代入式(1)中,可以得到子域E0中的位移試函數(shù)uh與目標(biāo)子域的支持域節(jié)點(diǎn)群自由度uE之間的關(guān)系如下

矩陣N為子域E0內(nèi)任意一點(diǎn)x=[x,y]T處的位移形函數(shù),其表達(dá)式如下

1.2 碎點(diǎn)法的弱形式動(dòng)量方程

在求解域 Ω 中,強(qiáng)形式動(dòng)量方程和邊界條件為:

式中 ρ為密度,是加速度向量,是施加的體力向量,σij是應(yīng)力張量.邊界條件中的 Γt和Γu分別代表施加應(yīng)力邊界條件和位移邊界條件的外邊界,ni是邊界上的外法線單位向量,和分別為在外邊界Γt和Γu上施加的面力和位移.

以虛位移 δui作為檢驗(yàn)函數(shù),采用伽遼金法推導(dǎo)弱形式動(dòng)量方程如下

對式(8) 括號內(nèi)第一項(xiàng)使用分部積分和高斯公式可得

式中 ?E表示子域E的邊界,nj為邊界?E上單位外法向量.

在碎點(diǎn)法模型中,任意子域E都需滿足式(9),對所有子域上的方程進(jìn)行累加可得

定義 Γ 表示所有子域邊界的集合,則Γh=ΓΓt-Γu表示所有相鄰子域間內(nèi)邊界的集合.假定任意相鄰子域E1和E2的公共邊界為e,對于任意物理量w,其在邊界e上的跳躍算子[]和平均算子{ }定義為

將跳躍算子和平均算子以及面力邊界條件引入式(10)的第四項(xiàng)中可得

式中 ηh內(nèi)邊界 Γh處的懲罰系數(shù),用于施加相鄰子域間的弱形式位移連續(xù)性約束,ηu為外邊界 Γu處的懲罰系數(shù),用于施加弱形式位移邊界條件,lh和lu分別為 Γh和Γu處子域界面的特征長度.

合并式(14)中等號右邊的第三和第四項(xiàng),并引入內(nèi)部界面數(shù)值通量的定義如下

最終推導(dǎo)得到包含數(shù)值通量的碎點(diǎn)法弱形式動(dòng)量方程如下

式中等號右邊第三項(xiàng)成為數(shù)值通量修正項(xiàng),可以看出在碎點(diǎn)法弱形式方程中,相鄰子域之間的相互作用僅由數(shù)值通量修正項(xiàng)決定,因此式(15)中定義的數(shù)值通量包含相鄰子域間相互作用面力的物理含義.

1.3 碎點(diǎn)法弱形式方程的離散

1.3.1 離散格式的弱形式方程離散格式的弱形式方程

為了便于表達(dá),以下推導(dǎo)過程均采用Voigt 標(biāo)記法的矩陣表達(dá)式.

由位移試函數(shù)與支持域節(jié)點(diǎn)自由度之間的關(guān)系式(7)可得速度和加速度的試探函數(shù)表達(dá)式如下

由于試探解與檢驗(yàn)函數(shù)采用相同形函數(shù),對虛位移 δu有

為子域E0的應(yīng)變形函數(shù).

將式(7)及式(17)～式(20)代入式(16)中,對碎點(diǎn)法求解域中所有子域進(jìn)行組裝,并在等號兩邊同時(shí)消去全局虛位移向量 δuE,可得到離散形式的碎點(diǎn)法動(dòng)量方程

其中,矩陣M是全局質(zhì)量矩陣,是全局加速度矩陣,fext是全局等效節(jié)點(diǎn)外力矩陣,fint是全局等效節(jié)點(diǎn)內(nèi)力矩陣,其表達(dá)式如下

式中,為子域的體力矩陣,為基于應(yīng)力邊界條件施加的面力,為子域的應(yīng)力,是內(nèi)邊界上的數(shù)值通量,是基于內(nèi)邊界單位法向量將子域應(yīng)力轉(zhuǎn)換成內(nèi)邊界面力的投影矩陣,為子域節(jié)點(diǎn)處的位移,為基于位移邊界條件施加的位移.

1.3.2 引入裂紋對弱形式方程的處理

碎點(diǎn)法使用的是分片連續(xù)的弱形式方程,因此易于在任意內(nèi)部界面處顯式引入裂紋.在離散格式的碎點(diǎn)法弱形式方程式(22)中,任意相鄰子域之間的聯(lián)系包含兩個(gè)部分,一部分是數(shù)值通量上的聯(lián)系,如上所述,數(shù)值通量t*具有相鄰子域相互所用面力的物理含義;另一部分是位移形函數(shù)上的聯(lián)系,因?yàn)樗辄c(diǎn)法的形函數(shù)是基于支持域節(jié)點(diǎn)群計(jì)算的,而相鄰子域互為對方的支持域.

由于數(shù)值通量的物理含義是相鄰子域的相互作用面力,可以自然地基于數(shù)值通量定義內(nèi)界面斷裂準(zhǔn)則.本文使用最大拉應(yīng)力準(zhǔn)則作為斷裂準(zhǔn)則,即當(dāng)碎點(diǎn)法模型中任意內(nèi)界面處的法向數(shù)值通量達(dá)到材料的拉伸強(qiáng)度時(shí),則判斷該內(nèi)界面發(fā)生斷裂.當(dāng)然,本文僅以最大拉應(yīng)力準(zhǔn)則為例,其他斷裂準(zhǔn)則也可以用作碎點(diǎn)法的內(nèi)界面斷裂判據(jù).

對于碎點(diǎn)法模型中的任意內(nèi)部邊界,當(dāng)邊界的數(shù)值通量滿足斷裂準(zhǔn)則時(shí),就在該邊界處引入裂紋.如上所述,任意相鄰子域之間的聯(lián)系包括兩個(gè)部分,因此需分別進(jìn)行修正.首先,對于斷裂的內(nèi)邊界,移除弱形式方程中該邊界處的數(shù)值通量修正項(xiàng),即將數(shù)值通量置零t*=0,就相當(dāng)于消除了數(shù)值通量部分的相互聯(lián)系,如圖3(b)所示.然后,對該內(nèi)界面兩端子域各自的支持域進(jìn)行修正,分別從各自的支持域中將對方移除,就相當(dāng)于消除了形函數(shù)部分的相互聯(lián)系.

圖3 碎點(diǎn)法裂紋的引入Fig.3 Steps of introducing crack in FPM

通過以上步驟,就可以在碎點(diǎn)法模型中任意內(nèi)部邊界處引入顯式裂紋.值得一提的是,采用上述方法在碎點(diǎn)法模型中引入裂紋時(shí),只修正了裂紋兩端子域各自的支持域和移除了該內(nèi)界面處的數(shù)值通量,這些修正只影響裂紋附近局部區(qū)域,對整體弱形式方程的影響較小,且不需要進(jìn)行網(wǎng)格重構(gòu),也不會額外增加整體模型的自由度.由此可見,碎點(diǎn)法只需要通過簡單的修正就可以在模型中的任意內(nèi)邊界處顯式引入裂紋,該操作僅對弱形式方程造成局部影響,且不增加自由度,因此碎點(diǎn)法可以簡單、高效地模擬結(jié)構(gòu)中任意位置的裂紋萌生和任意形貌的裂紋擴(kuò)展,在模擬斷裂、破碎等極端問題時(shí)具有一定的優(yōu)勢.

1.4 碎點(diǎn)法的顯式動(dòng)力學(xué)求解

本文采用中心差分法來對碎點(diǎn)法弱形式動(dòng)量方程進(jìn)行求解,中心差分法的時(shí)間軸如圖4 所示.

圖4 中心差分法時(shí)間軸示意圖Fig.4 Time axis of the Verlet method

將式(24)和式(25)整理可得

這種形式的中心差分法也通常被稱作蛙跳格式.基于該方法,定義動(dòng)力學(xué)碎點(diǎn)法的求解過程如下:

(1) 由式(23)計(jì)算tn時(shí)刻的質(zhì)量和節(jié)點(diǎn)力矩陣

(7) 重復(fù)上述(1)～(6)步,直至最后一個(gè)時(shí)間步.

由于中心差分法是條件穩(wěn)定的,因此為了保證計(jì)算的穩(wěn)定性,每一個(gè)時(shí)間步的時(shí)間增量Δt必須小于臨界時(shí)間步長Δtcr.

下面討論中心差分法的穩(wěn)定性問題.為了維持算法的穩(wěn)定性,時(shí)間步長不得超過臨界時(shí)間步長Δtcr,一般情況下取

式中 α 是安全系數(shù),碎點(diǎn)法參考有限元的安全系數(shù)取值,取0.8 ≤α ≤0.98,臨界時(shí)間步長的計(jì)算公式如下

式中l(wèi)e代表碎點(diǎn)法模型中子域的特征長度,c是線彈性材料的絕熱聲速.

由于碎點(diǎn)法的弱形式方程引入了數(shù)值通量修正,用于保證方法的一致性,而數(shù)值通量修正項(xiàng)中包含懲罰系數(shù) ηh,如式(15)所示.數(shù)值通量會使得整個(gè)系統(tǒng)產(chǎn)生一定的數(shù)值振蕩,而數(shù)值通量產(chǎn)生的作用力的大小取決于懲罰系數(shù) ηh,因而懲罰系數(shù)會對臨界時(shí)間步長產(chǎn)生影響,較小的時(shí)間步長可以逐漸消除數(shù)值通量帶來的振蕩的影響.

對碎點(diǎn)法中中心差分法的臨界時(shí)間步長計(jì)算公式進(jìn)行修正得到

使用上式對碎點(diǎn)法模型中的每一個(gè)子域進(jìn)行計(jì)算,使用最小的臨界時(shí)間步長作為整體模型的臨界時(shí)間步長.

2 碎點(diǎn)法動(dòng)力學(xué)算例

2.1 拉力作用下的平板應(yīng)力波傳播

為了驗(yàn)證動(dòng)力學(xué)碎點(diǎn)法模擬應(yīng)力波傳播的能力,首先應(yīng)用動(dòng)力學(xué)碎點(diǎn)法程序求解如圖5 所示平板的應(yīng)力波傳播問題.參考文獻(xiàn)中的設(shè)置[40],該算例中使用矩形薄板作為研究對象,矩形板的面內(nèi)尺寸為長L=4 m,寬D=2 m.平板材料屬性為,楊氏模量E=80kPa,泊松比 μ=0,而材料密度 ρ=1 kg/m3.邊界條件如圖5 所示,矩形板左端固支,右端受均布拉伸載荷P(t),載荷大小隨試件的變化關(guān)系如圖6所示,仿真的總時(shí)間設(shè)置為0.3 s.由于分析對象為薄板,因此可將該問題等效為平面應(yīng)力問題,碎點(diǎn)法模型中使用40×20個(gè)均勻分布的節(jié)點(diǎn)對矩形板進(jìn)行離散.由于該問題中材料的泊松比設(shè)置為零,因此本質(zhì)上這是一個(gè)一維應(yīng)力波傳播問題.

圖5 拉力作用下的矩形平板模型Fig.5 The model for rectangular plate under tension

圖6 載荷 P(t) 隨時(shí)間的關(guān)系Fig.6 The curve for the load-time correlation

對于此問題,文獻(xiàn)[40]給出了自由端A點(diǎn)的水平位移uA、中點(diǎn)B點(diǎn)的水平位移位移uB、中點(diǎn)B點(diǎn)的水平正應(yīng)力、以及固定端C點(diǎn)的水平正應(yīng)力的精確解.本文以上述數(shù)據(jù)為參考對碎點(diǎn)法的計(jì)算結(jié)果進(jìn)行驗(yàn)證,同時(shí)與有限元結(jié)果進(jìn)行對照.如圖7 所示,在碎點(diǎn)法模型中的相同位置讀取位移和應(yīng)力信息,峰值處的相對誤差如表1和表2 所示.碎點(diǎn)法的仿真結(jié)果與有限元及理論解吻合良好,從而驗(yàn)證了本文所提出的動(dòng)力學(xué)碎點(diǎn)法在模擬一維應(yīng)力波傳播問題方面的有效性.

圖7 拉力作用下平板應(yīng)力波傳播問題精確解及FPM 結(jié)果Fig.7 Exact solutions and FPM’s results of the displacement and stress at different points

圖7 拉力作用下平板應(yīng)力波傳播問題精確解及FPM 結(jié)果(續(xù))Fig.7 Exact solutions and FPM’s results of the displacement and stress at different points(continued)

表1 有限元、碎點(diǎn)法的位移峰值與理論解的比較Table 1 Relative error of maximum displacement obtained from FEM and FPM

表2 有限元、碎點(diǎn)法的應(yīng)力峰值與理論解的比較Table 2 Relative error of maximum stress obtained from FEM and FPM

2.2 自由端剪切載荷作用下懸臂梁動(dòng)態(tài)響應(yīng)

接著使用動(dòng)力學(xué)碎點(diǎn)法程序?qū)ψ杂啥思羟休d荷作用下的懸臂梁動(dòng)態(tài)響應(yīng)進(jìn)行仿真,通過這個(gè)算例來驗(yàn)證所提出的方法在模擬二維應(yīng)力波傳播問題方面的能力.該算例中使用圖8(a)所示的懸臂梁作為研究對象,梁的尺寸為,長度L=48 m,高度H=12 m.材料的彈性模量E=30MPa,泊松比為0.3,密度ρ=1.0kg/m3.邊界條件設(shè)置如圖8(a)所示,懸臂梁左端面固支,右端面受向下的瞬態(tài)剪切載荷P=100Pa,載荷作用時(shí)間為0～0.5 s,總仿真時(shí)間為 2 s.本算例同樣使用平面應(yīng)力假設(shè),碎點(diǎn)法模型使用 4 8×12 個(gè)均勻分布的節(jié)點(diǎn)對懸臂梁進(jìn)行離散,如圖8(b).從模型中隨機(jī)選取一個(gè)節(jié)點(diǎn)A點(diǎn),用于讀取位移和應(yīng)力仿真結(jié)果.本算例中所設(shè)置的幾何、材料參數(shù)僅用于算法驗(yàn)證,不代表任何具體的物理對象.

圖8 (a)懸臂梁示意圖和(b)計(jì)算模型節(jié)點(diǎn)分布Fig.8 (a) Configuration of the beam and(b) the scatter of subdomains

同時(shí),還在商業(yè)有限元軟件ABAQUS 中建立了這個(gè)算例的有限元模型,用有限元模擬結(jié)果對碎點(diǎn)法程序進(jìn)行驗(yàn)證.有限元模型的網(wǎng)格與碎點(diǎn)法模型的子域完全一致.

在圖9～ 圖11 中分別展示了懸臂梁加載端(右端)中點(diǎn)的豎直位移、固定端(左端)中點(diǎn)各應(yīng)力分量和A點(diǎn)水平正應(yīng)力的有限元方法FEM和碎點(diǎn)法計(jì)算結(jié)果.圖11 中僅繪制A點(diǎn)的水平正應(yīng)力曲線是因?yàn)樵擖c(diǎn)處的其他應(yīng)力分量基本為零.可以看出,碎點(diǎn)法預(yù)測的應(yīng)力和位移狀態(tài)與有限元結(jié)果吻合良好,兩種方法得到的應(yīng)力、位移響應(yīng)頻率特征一致.各圖中,兩種方法得到的位移曲線幾乎完全重合,應(yīng)力曲線峰值的相對誤差在 2.5% 以內(nèi),其余部分幾乎完全重合.由此可見,本文所提出的動(dòng)力學(xué)碎點(diǎn)法程序可以有效模擬二維應(yīng)力波傳播問題.

圖9 自由端中點(diǎn)豎直方向位移響應(yīng)曲線Fig.9 Vertical displacement curve of the mid-point at the free end of the beam

圖10 固定端中點(diǎn)應(yīng)力響應(yīng)曲線Fig.10 Curves of stress components of the mid-point at the fixed end of the beam

圖11 懸臂梁A 點(diǎn)處應(yīng)力瞬態(tài)響應(yīng)曲線Fig.11 Stress-time curve of point A

2.3 含裂紋金屬板低速?zèng)_擊破壞行為仿真

最后,使用動(dòng)力學(xué)碎點(diǎn)法程序?qū)σ粋€(gè)經(jīng)典動(dòng)態(tài)斷裂問題進(jìn)行模擬仿真,以驗(yàn)證該方法模擬裂紋動(dòng)態(tài)擴(kuò)展的能力.該算例基于 Kalthoff等[41-42]在1987 年和2000年所做的一系列試驗(yàn)工作,試驗(yàn)試件如圖12(a)所示,為長和寬分別為100mm和200mm的矩形金屬板,板上包含著兩條對稱的水平初始裂紋,裂紋間距為50mm,初始裂紋長度為50mm.試驗(yàn)使用一塊底面半徑為25 mm 的圓柱形子彈以33 m/s的速度沖擊兩條裂紋之間的區(qū)域[42].試驗(yàn)觀察到含裂紋金屬板在子彈沖擊下的主要失效模式為脆性斷裂,裂紋擴(kuò)展路徑與初始裂紋的夾角約為70°,如圖12(b)所示.

圖12 (a)含裂紋金屬板沖擊試驗(yàn)裝置示意圖和(b)試驗(yàn)觀測到的裂紋擴(kuò)展路徑[43]Fig.12 (a) Experimental setup for the Kalthoff plate impact test and(b) the experimentally-observed crack path[43]

根據(jù)該問題的對稱性,建立二分之一模型以減少計(jì)算量,計(jì)算模型如圖13(a)所示,圖中紅色線條代表初始裂紋.文獻(xiàn)[33]研究了碎點(diǎn)法節(jié)點(diǎn)分布方式給計(jì)算結(jié)果帶來的影響,由于數(shù)值通量修正的引入,碎點(diǎn)法使用均勻分布節(jié)點(diǎn)或隨機(jī)分布節(jié)點(diǎn)均可以得到準(zhǔn)確的結(jié)果,驗(yàn)證了碎點(diǎn)法的魯棒性.基于該特點(diǎn)及試驗(yàn)中觀測到的裂紋主要擴(kuò)展區(qū)域,在布置節(jié)點(diǎn)時(shí),矩形板右上區(qū)域布置的子域節(jié)點(diǎn)更加密集,在其他區(qū)域子域節(jié)點(diǎn)稍微稀疏.這種方式一方面可以更精細(xì)地捕捉裂紋的擴(kuò)展路徑,另一方面可以節(jié)省計(jì)算資源.本文的工作中,鄰近子域節(jié)點(diǎn)群的選取為所有與該子域共邊界的子域節(jié)點(diǎn),并未因節(jié)點(diǎn)分布密度不同而進(jìn)行調(diào)整.

試驗(yàn)中金屬板由型號18 Ni1900的鋼材制成,碎點(diǎn)法仿真中采用文獻(xiàn)[44]中的材料參數(shù),楊氏模量E=190GPa,泊松比 μ=0.3,密度 ρ=8 t/m3.根據(jù)文獻(xiàn)[45],取內(nèi)界面破壞的拉伸強(qiáng)度為1773 MPa.同時(shí),還建立了使用更加致密節(jié)點(diǎn)的碎點(diǎn)法模型,如圖13(b)所示,用于分析網(wǎng)格密度對仿真結(jié)果的影響,并驗(yàn)證圖13(a)所示的有效性.

圖13 計(jì)算模型:(a) 6516 個(gè)子域節(jié)點(diǎn)和(b) 11 995 個(gè)子域節(jié)點(diǎn)Fig.13 Domain discretization with(a) 6516 FPM subdomains and(b)11 995 FPM subdomains

由于沖擊載荷作用時(shí)間短,可認(rèn)為子彈速度在沖擊過程中幾乎沒有降低,因此可以在沖擊區(qū)域用速度邊界條件等效代替子彈的沖擊載荷.根據(jù)文獻(xiàn)[45-46],試驗(yàn)中所使用的子彈和試件由相同材料制備而成,二者擁有相同的彈性阻抗,因此加載速度取試驗(yàn)中子彈速度的一半 1 6.5 m/s.

使用碎點(diǎn)法仿真計(jì)算得到的靜水壓力云圖如圖14 所示,可以看到圖14(a)中,在 8 μs 時(shí),由于矩形板左側(cè)載荷的作用,應(yīng)力波先是沿著矩形板下方初始裂紋與對稱面之間的區(qū)域,而矩形板其他區(qū)域暫時(shí)未受到擾動(dòng),相應(yīng)的靜水壓力均為零.到12 μs時(shí),裂紋尖端處靜水壓力較大,矩形板其他位置的靜水壓力幾乎為零,如圖14(b)所示,這是因?yàn)槌跏剂鸭y的存在導(dǎo)致了應(yīng)力波逐漸在初始裂紋尖端處聚集,產(chǎn)生應(yīng)力集中的現(xiàn)象,仿真的結(jié)果與實(shí)際的物理機(jī)制吻合.

由于裂紋尖端應(yīng)力集中,裂紋尖端區(qū)域的內(nèi)邊界達(dá)到強(qiáng)度而開始萌生裂紋,并沿著與初始裂紋夾角約70°角的方向擴(kuò)展,如圖14(c)所示,此時(shí)應(yīng)力集中的位置也隨之發(fā)生變化,應(yīng)力集中現(xiàn)象發(fā)生在新的裂紋尖端處,而已生成裂紋的區(qū)域應(yīng)力逐漸減小,仿真的結(jié)果吻合物理機(jī)制.

圖14 靜水壓云圖Fig.14 Configuration with hydrostatic pressure

圖14 靜水壓云圖(續(xù))Fig.14 Configuration with hydrostatic pressure(continued)

不同節(jié)點(diǎn)分布的仿真的裂紋擴(kuò)展路徑在圖15中畫出.對于仿真結(jié)果1,在最初時(shí)刻,裂紋以大約70°角開始擴(kuò)展.隨著裂紋擴(kuò)展到水平位置60mm左右,夾角出現(xiàn)偏轉(zhuǎn),此時(shí)角度約為64°;在水平位置80mm 處左右,裂紋近似沿垂直方向擴(kuò)展,直至上邊緣附近,又重新以約70°角擴(kuò)展;從初始裂紋尖端到裂紋最終位置的平均角度約為66°.而對于仿真結(jié)果2,裂紋幾乎一致沿著約70°的角度擴(kuò)展,僅在接近板上邊緣附近出現(xiàn)了角度增大的情況.總體而言,兩種不同節(jié)點(diǎn)密度的模型預(yù)測的裂紋擴(kuò)展形貌都和試驗(yàn)觀測基本一致.其中模型2 使用更致密的節(jié)點(diǎn),因此仿真結(jié)果與試驗(yàn)結(jié)果吻合更加良好;模型1 的自由度顯著少于模型2,計(jì)算效率更高,且能夠給出關(guān)鍵性的應(yīng)力傳播和裂紋擴(kuò)展特性,也滿足動(dòng)態(tài)裂紋擴(kuò)展預(yù)測和分析的需求.

圖15 FPM 仿真裂紋擴(kuò)展路徑Fig.15 Crack propagation from FPM simulation

下面繪制裂紋擴(kuò)展的速度并與文獻(xiàn)[45]的仿真結(jié)果對比,討論碎點(diǎn)法仿真得到的裂紋擴(kuò)展速度的響應(yīng)特征.根據(jù)瑞利波速計(jì)算公式

式中cs為剪切波傳播速度,可以計(jì)算得到瑞利波速為cR=2799 m/s.

對裂紋長度隨時(shí)間的變化數(shù)據(jù)進(jìn)行差分求導(dǎo)處理可求得裂紋擴(kuò)展速度,如圖16 所示,裂紋擴(kuò)展速度在25 μs 左右時(shí)開始增長,之后趨于穩(wěn)定,在56 μs左右開始出現(xiàn)下降.碎點(diǎn)法得到的裂紋擴(kuò)展速度略低于文獻(xiàn)仿真結(jié)果,但總體趨勢與文獻(xiàn)[45]中數(shù)據(jù)吻合良好.

圖16 裂紋擴(kuò)展速度-時(shí)間曲線Fig.16 Curve of crack propagation speed versus time

以上算例表明,本文所提出的動(dòng)力學(xué)碎點(diǎn)法能夠模擬典型的動(dòng)態(tài)裂紋擴(kuò)展行為,仿真得到的裂紋形貌和裂紋擴(kuò)展速度與試驗(yàn)觀測結(jié)果吻合良好.考慮到碎點(diǎn)法可以在任意內(nèi)界面處判斷界面狀態(tài),當(dāng)界面破壞時(shí)只需要簡單的修正就可以顯式引入裂紋,且引入裂紋不會增加模型整體的自由度,因此動(dòng)力學(xué)碎點(diǎn)法具備模擬大規(guī)模復(fù)雜裂紋萌生和動(dòng)態(tài)擴(kuò)展問題的能力.

3 結(jié)論及展望

本文根據(jù)碎點(diǎn)法的基本思想和原理,推導(dǎo)了碎點(diǎn)法弱形式動(dòng)量方程、建立了顯式動(dòng)力學(xué)碎點(diǎn)法求解格式并編制了相應(yīng)的計(jì)算程序,還使用經(jīng)典算例驗(yàn)證了動(dòng)力學(xué)碎點(diǎn)法模擬應(yīng)力波傳播和動(dòng)態(tài)裂紋擴(kuò)展問題的能力.

碎點(diǎn)法具有如下特色:(1)碎點(diǎn)法是一種基于伽遼金弱形式的無網(wǎng)格法,具有較好的算法穩(wěn)定性;(2)采用基于節(jié)點(diǎn)的位移形函數(shù),因此具有抵抗子域畸變的能力;(3)碎點(diǎn)法的形函數(shù)是分片連續(xù)的,因此可以使用多項(xiàng)式形式的形函數(shù),并使用簡單的高斯積分進(jìn)行求解;(4)在弱形式方程中引入數(shù)值通量修正項(xiàng),從而保證放棄了內(nèi)界面位移連續(xù)性的碎點(diǎn)法的一致性和穩(wěn)定性;(5)內(nèi)界面數(shù)值通量具有相鄰子域間相互作用面力的物理含義,因此可用于定義內(nèi)界面的斷裂準(zhǔn)則;(6)只需要簡單的操作就可以在碎點(diǎn)法模型中顯式引入裂紋,即移除對應(yīng)內(nèi)界面的數(shù)值通量修正,并修改界面兩端子域的支持域,從而分別相鄰子域間數(shù)值通量和形函數(shù)兩方面的相互聯(lián)系;(7)子域節(jié)點(diǎn)具有質(zhì)量、體積、速度等明確的物理屬性,便于顯式動(dòng)力學(xué)計(jì)算.

隨后通過三個(gè)典型算例驗(yàn)證了動(dòng)力學(xué)碎點(diǎn)法對應(yīng)力波傳播及動(dòng)態(tài)裂紋預(yù)測的預(yù)測能力.前兩個(gè)算例驗(yàn)證了碎點(diǎn)法模擬拉伸和剪切應(yīng)力波傳播的有效性和準(zhǔn)確性,所得結(jié)果與精確解和有限元結(jié)果吻合良好.第三個(gè)算例采用經(jīng)典動(dòng)態(tài)斷裂算例來驗(yàn)證動(dòng)力學(xué)碎點(diǎn)法模擬裂紋動(dòng)態(tài)擴(kuò)展過程的能力,計(jì)算結(jié)果與文獻(xiàn)中的試驗(yàn)結(jié)果吻合良好,體現(xiàn)了碎點(diǎn)法模擬動(dòng)態(tài)斷裂的問題的有效性.本文所提出的動(dòng)力學(xué)碎點(diǎn)法理論及程序?yàn)閯?dòng)態(tài)斷裂問題的研究提供了一種簡單、高效的數(shù)值仿真方法.

最后需要指出的是,本文并未涉及如裂紋分叉、交匯等問題,但碎點(diǎn)法易于顯式引入裂紋的特點(diǎn)使得其在這類問題的仿真中具有很大潛力.這些問題背后的復(fù)雜物理機(jī)制意味著需要做進(jìn)一步工作,如開發(fā)更加細(xì)化的界面模型,確定合適的斷裂準(zhǔn)則,以及建立這些準(zhǔn)則與界面數(shù)值通量的聯(lián)系等,這也將是作者團(tuán)隊(duì)下一步的研究方向.