等幾何分析與近場動力學(xué)耦合的殼結(jié)構(gòu)裂紋擴展仿真算法

夏陽， 王洪帥， 鄭國君， 申國哲

(大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室 汽車工程學(xué)院，遼寧 大連 116024)

0 引 言

斷裂仿真對殼結(jié)構(gòu)的安全評估具有重要意義。近年來，斷裂力學(xué)的數(shù)值仿真取得快速發(fā)展，為業(yè)界提供多種可行的解決方案，如內(nèi)聚力單元法、擴展有限元方法、擴展等幾何方法、無網(wǎng)格法、相場方法、分子動力學(xué)方法和近場動力學(xué)方法等。目前，斷裂仿真存在的問題是如何對結(jié)構(gòu)中連續(xù)場和非連續(xù)場進(jìn)行統(tǒng)一描述，并且在計算中保證計算效率。為解決這一問題，斷裂力學(xué)非連續(xù)模型與連續(xù)模型相結(jié)合是一個值得研究的方向。

等幾何分析(isogeometric analysis, IGA)直接將精確幾何用于數(shù)值分析，連通幾何設(shè)計與數(shù)值分析之間的鴻溝，適用于殼結(jié)構(gòu)分析的等幾何列式不斷被提出，包括無轉(zhuǎn)動的Kirchhoff-Love(K-L)殼模型、帶轉(zhuǎn)動自由度的Reissner-Mindlin(R-M)殼模型、實體殼等?；诹鸭y擴展問題，將擴展有限元中用到的單位分解法引入IGA中，得到基于K-L理論和R-M理論的擴展IGA(extended IGA, XIGA)方法。XIGA能夠模擬裂紋附近的應(yīng)力場，但是需要額外的準(zhǔn)則處理裂紋成核和擴展，在處理裂紋分叉和合并方面存在一定困難。另外一種處理斷裂問題的耦合方法是IGA與相場方法耦合。相場方法是建立在Griffith能量斷裂準(zhǔn)則之上的，該準(zhǔn)則是定量描述斷裂擴展的經(jīng)典準(zhǔn)則，其核心是假設(shè)裂紋擴展過程中的總能量(表面能和體積能)保持最小。AMBATI等采用等幾何實體殼單元與相場方法耦合的方式處理殼結(jié)構(gòu)的脆性和韌性斷裂。KIKIS等將考慮剪切變形的等幾何R-M殼與相場方法相結(jié)合，避免體積離散化。為能夠得到合適的結(jié)果，相場方法必須在裂紋擴展處進(jìn)行網(wǎng)格細(xì)化以滿足收斂的要求，導(dǎo)致其需要消耗大量的計算時間。

近場動力學(xué)(peridynamics, PD)理論由SILLING提出，用于處理非連續(xù)場問題。學(xué)者們提出很多用于板殼結(jié)構(gòu)問題的PD模型。TAYLOR等通過將3D模型降為2D，提出PD薄板模型，該模型能夠捕獲面外彎曲變形，但是只能處理具有特定泊松比的材料。DIYAROGLU等提出考慮剪切變形的薄/厚梁和板的PD列式，其運動方程可以簡化為經(jīng)典Timoshenko梁和Mindlin板方程。SHEN等基于微梁鍵提出6自由度Euler梁和K-L殼模型，通過插值方法獲得鍵的變形、建立鍵的微勢能，從而得到鍵的微模量，該模型不限制泊松比。以此為基礎(chǔ)，考慮剪切變形，進(jìn)一步發(fā)展出PD R-M殼模型。

由于PD是一種非局部理論，每一個點與其近場域內(nèi)的一群點相互作用，計算效率比經(jīng)典的有限元或IGA更低。采用非局部理論和經(jīng)典連續(xù)介質(zhì)力學(xué)理論耦合的方式，其計算效率有望得到提高。將IGA與PD耦合是進(jìn)行裂紋擴展仿真的一種新的重要方法。MADENCI等提出IGA與PD耦合處理平面問題的方法，將PD點加入到等幾何模型中，PD點的自由度與等幾何節(jié)點相關(guān)聯(lián)，但是PD點不出現(xiàn)在最終的方程中。XIA等提出用于平面彈性問題的IGA-PD耦合模型，使用等幾何控制點作為PD點，根據(jù)力平衡方法將IGA方程與PD方程相疊加得到待求解線性方程，是一種簡單且高效的求解方法。此外，鄭國君等使用神經(jīng)網(wǎng)絡(luò)預(yù)測裂紋擴展路徑，快速且準(zhǔn)確。

殼的光滑曲面特征使其適于采用IGA方法進(jìn)行模擬。為對殼結(jié)構(gòu)進(jìn)行裂紋擴展分析，本文提出一種用于板殼結(jié)構(gòu)斷裂分析的耦合模型IGA-PD。采用IGA模型保證幾何精確性，將PD模型用于裂紋擴展區(qū)域。為直觀顯示裂紋擴展過程，采用非連續(xù)伽遼金(discontinuous Galerkin, DG)方程，提出DG有限元PD板殼單元列式。與粒子點法PD列式相比，采用DG有限元PD板殼單元可以使用高斯點積分，提高計算精度。在耦合過程中，采用力/力矩平衡原理得到耦合模型。將IGA作為邊界可以消除PD的表面效應(yīng)問題，提高PD模型性能。

1 R-M殼的IGA和微梁鍵模型

1.1 R-M殼IGA模型

采用由三維彈性模型退化得到的R-M殼模型，其幾何模型示意見圖1。殼變形前采用直法線假設(shè)，其物理場采用非均勻有理B樣條(non-uniform rational B-spline, NURBS)插值，表達(dá)式為

圖 1 R-M殼模型示意

(,,)=(,)+·(,)=

(1)

式中：和為節(jié)點參數(shù),可分別為由節(jié)點向量=[…++1]和=[…++1]定義；∈[-1,1]為殼的厚度；=[]，即殼上控制點的位置矢量；(,)為與控制點有關(guān)的二維張量積NURBS基函數(shù)；為控制點處殼體中面的法向量，=[，1，2,3]，采用Greville橫標(biāo)點計算。

基于等參概念，采用與幾何場相同的插值函數(shù)，殼的位移表示為

(,,)=

(2)

對于階基函數(shù)，完全積分需要+1個高斯點。位移場與轉(zhuǎn)動場不一致會導(dǎo)致剪切閉鎖，因此使用個高斯點的減縮積分方案避免這一問題。

通過最小位能原理，最終得到等幾何R-M殼單元的剛度矩陣

=?d=

(3)

在圖1中給定一點,其坐標(biāo)為(,,)，根據(jù)式(2)，點的速度可以通過點計算得到，即

=+

(4)

其中：

(5)

點處體積微元的動能為

(6)

根據(jù)直法線假設(shè)，點相對于點的位置表示為=[0 0]。將式(4)代入式(6)，得到采用中面表示的動能方程，即

(7)

式中：為材料密度；

(8)

(9)

采用等幾何形函數(shù)進(jìn)行離散，對式(7)在單元域上進(jìn)行積分，得到一致質(zhì)量矩陣，即

(10)

1.2 適用于R-M殼的微梁鍵

近場動力學(xué)是一種使用積分方程求解連續(xù)和非連續(xù)場的理論。該理論假設(shè)粒子與周圍一定距離范圍內(nèi)的其他粒子存在相互作用，損傷發(fā)生在粒子水平上。參考構(gòu)型中2個粒子的相對位置為

=′-

(11)

相對位移為

=(′,)-(,)

(12)

這2個相互作用的粒子之間形成PD鍵。粒子作用在粒子上的力可由微勢能求導(dǎo)得到，即

(13)

微勢能是單根鍵上的能量，其單位是單位體積的平方，所以給定點的單位體積的能量即局部應(yīng)變能密度

(14)

系數(shù)12是指每根鍵的每個端點上只有鍵上12的能量?；诮鼒鰟恿W(xué)理論的微梁鍵，將PD鍵視為類似于經(jīng)典Timoshenko梁模型，承受軸向變形、扭轉(zhuǎn)變形、彎曲變形和橫向剪切變形，稱之為微梁鍵，見圖2。

圖 2 PD微梁鍵示意

每個微梁鍵的端點有3個平動自由度(,,)和3個轉(zhuǎn)動自由度(,,)，通過插值方式可獲得鍵的變形。鍵在局部坐標(biāo)系下的位移向量為

(15)

根據(jù)文獻(xiàn)[14-15]，為描述殼結(jié)構(gòu)的變形，微梁鍵上一點的變形可以分成面內(nèi)變形和橫向變形兩部分，即

,Shell=,PlaneShell+,TranShell

(16)

面內(nèi)變形微勢能由端點平動產(chǎn)生的軸向變形能和端點轉(zhuǎn)動產(chǎn)生的橫向變形能組成，則

(17)

式中：a,為鍵在拉伸和壓縮作用下的微模量；b,為與面內(nèi)彎曲相關(guān)的微模量；為鍵的長度。

由殼橫向變形相關(guān)的微勢能可得

(18)

式中：為與轉(zhuǎn)動相關(guān)的微模量；b,和s,分別為由于鍵繞軸轉(zhuǎn)動產(chǎn)生的與彎曲和橫向剪切相關(guān)的微模量。

根據(jù)應(yīng)變能密度守恒原理，微梁鍵基殼模型的應(yīng)變能密度必須與彈性殼理論的應(yīng)變能密度相等。通過插值離散微梁鍵模型，得到微模量的表達(dá)式為

(19)

最終，局部坐標(biāo)系下微梁鍵的力密度與位移的關(guān)系表示為

=

(20)

(21)

(22)

的具體細(xì)節(jié)可參閱文獻(xiàn)[14-15]。

2 非連續(xù)伽遼金PD殼

在經(jīng)典PD模型中，物理域被離散成粒子，稱為無網(wǎng)格PD離散方法；近幾年發(fā)展的另一種方法是將模型域離散成單元，再應(yīng)用DG方程。例如，REN等使用DG模型離散裂紋存在的非局部模型，使用有限元法離散連續(xù)模型域。與粒子法相比，DG方法可以在計算域內(nèi)使用更高精度的積分，可以更方便地計算應(yīng)力和應(yīng)變場。采用微梁鍵基PD殼模型計算的系統(tǒng)總能量為

(23)

(24)

式中：為全局坐標(biāo)系下點和的位移；為全局坐標(biāo)系到鍵局部坐標(biāo)系下的旋轉(zhuǎn)矩陣。

使用DG方法時，需要先將模型離散成單元，本文采用由控制網(wǎng)形成的4節(jié)點單元，見圖3。為所有單元的共享節(jié)點插入新的節(jié)點，解除單元之間的節(jié)點關(guān)聯(lián)。這種單元的拆分過程可保證各個單元的獨立性，使裂紋擴展分析可以順利進(jìn)行。

圖 3 DG離散中的單元拆分示意

在PD殼的DG法中，單元是計算模型的基本組成部分。在列式中，單元的高斯點被視為PD節(jié)點，見圖4。

圖 4 DG單元的高斯點之間形成PD鍵

單元的高斯點被視為附著在單元上的PD節(jié)點。兩單元的任意2個高斯點之間形成鍵，如圖中鍵連接點和點，點的近場域由單元決定。單元中的4個PD點都采用單元形心形成的近場域，單元和單元之間的關(guān)聯(lián)通過鍵的相互作用建立，即

(25)

令和分別表示單元坐標(biāo)系與全局坐標(biāo)系之間位移向量的轉(zhuǎn)換矩陣，則

(26)

將式(24)～(26)代入式(23)，得

(27)

采用DG法的單元與單元之間的單剛矩陣可寫為

=

(28)

對所有單元的剛度矩陣進(jìn)行累加，即得到結(jié)構(gòu)的整體剛度矩陣。

3 基于力/力矩平衡原理的耦合模型

DG列式方法增加單元節(jié)點的數(shù)量，因此降低計算效率。為提高裂紋擴展仿真的計算效率，根據(jù)前文獲得的IGA殼模型和DG殼的PD模型，采用非局部模型和經(jīng)典連續(xù)介質(zhì)力學(xué)模型耦合的方法，得到殼結(jié)構(gòu)的動力學(xué)平衡方程，即

(29)

方程的每一行表示一個自由度所處的一種力力矩平衡狀態(tài)。處于一種邊界條件下的結(jié)構(gòu)，其所承受的等效載荷向量相等，即滿足=。當(dāng)采用IGA模型或PD模型求解時，2種模型得到的位移解應(yīng)該相同，即滿足=。因此，方程又可表示為

(30)

根據(jù)式(30)，可提出基于力力矩平衡的耦合求解方法。將計算域按等幾何控制點形成的網(wǎng)格劃分為IGA域與PD域，見圖5。

圖 5 耦合模型示意

由此，耦合剛度矩陣可表示為

(-())+()

(31)

式中：為對角矩陣，

(32)

且滿足

(33)

耦合后的動力學(xué)平衡方程表示為

(34)

由式(30)可以清楚地看出，式(34)中的每個方程都成立。節(jié)點各自由度對應(yīng)的力和力矩狀態(tài)是平衡的，因此這種方法稱為力力矩平衡耦合方法。

對于線彈性靜力學(xué)問題，可以通過求解線性方程獲得模型的位移解，即

=

(35)

對于裂紋動態(tài)擴展問題，可采用中心差分法求解。系統(tǒng)待求解的運動方程為

(36)

4 算 例

4.1 帶初始裂紋板的裂紋擴展

本算例被廣泛用于評估裂紋擴展仿真算法的準(zhǔn)確性。初始裂紋板模型及PD區(qū)域設(shè)置見圖6，板的材料參數(shù)為=72 GPa，=033，=2 440 kg/m，=135 J/m；作用在模型上下邊界的均布載荷為=12 MPa；將近場域大小設(shè)置為=3Δ。

圖 6 初始裂紋板模型,mm

計算得到不同時間步下的裂紋擴展路徑，見圖7。IGA-PD耦合模型模擬殼裂紋擴展是有效的。同時，根據(jù)裂紋擴展區(qū)域的分布，可選擇合適的區(qū)域作為PD區(qū)。對于裂紋擴展區(qū)域分布集中的問題，可以采用IGA-PD耦合算法提高裂紋擴展分析的計算效率。

圖 7 初始裂紋板的裂紋擴展路徑

4.2 平行預(yù)置裂紋板的裂紋擴展

為檢驗耦合算法在復(fù)雜裂紋條件下對裂紋擴展和合并的仿真能力，采用受位移載荷的預(yù)置平行裂紋板算例，見圖8?？紤]裂紋角度=0和=45°這2種情況，幾何參數(shù)=10 mm，=10 mm。模型參數(shù)=30 GPa，=033，=2 700 kg/m，=30 J/m；將近場域大小設(shè)置為=3Δ。模型上下邊界給定0.04 mm位移，分2 000步線性加載，時間步長Δ=5×10s，則位移加載速度為0.4 m/s。PD區(qū)域根據(jù)預(yù)制裂紋分布進(jìn)行布置。對于圖8中的傾斜裂紋，可設(shè)置包圍裂紋的區(qū)域為PD區(qū)域。裂紋擴展過程見圖9和10。2種角度的裂紋均在開始階段各自生長，內(nèi)側(cè)相互靠近并逐漸合二為一，外側(cè)向兩側(cè)擴展至邊緣。裂紋擴展路徑與文獻(xiàn)[23]中的結(jié)果一致。

圖 8 平行預(yù)置裂紋板模型,mm

圖 9 α=0°平行裂紋板的裂紋擴展路徑

圖 10 α=45°平行裂紋板的裂紋擴展路徑

采用=0的模型研究不同PD區(qū)域大小對計算效率的影響。所用測試設(shè)備參數(shù)為CPU Intel i5-7400 3.00 GHz，RAM 32 GB 2 400 MHz，硬盤TOSHIBA DT01ACA100 7200 HDD 1 TB。以劃定的PD區(qū)域面積占整個模型總面積的比例為橫坐標(biāo)，以采用耦合模型求解時間與采用純PD模型計算時間的比值為縱坐標(biāo)，不同PD區(qū)域占比的計算時間與純PD區(qū)域的計算時間的相對用時對比見圖11。由此可知，隨著PD區(qū)域所占比例的增加，計算時間也增加。采用20%區(qū)域為PD區(qū)域，其余區(qū)域為IGA模型區(qū)域時，所需計算時間為純PD模型的25%，說明與PD方法相比，使用IGA-PD耦合方法能夠提高計算效率。

圖 11 不同PD區(qū)域時的計算時間對比

4.3 受內(nèi)壓的圓柱殼

為檢驗耦合算法在一般殼結(jié)構(gòu)斷裂仿真方面的性能，采用圓柱作為研究對象，模擬其在受內(nèi)部壓強作用下的斷裂過程。帶預(yù)制裂紋的受內(nèi)壓圓柱殼模型見圖12。材料參數(shù)=70 GPa，=03，=1.5 J/m；模型承受的內(nèi)部壓強=1 MPa；設(shè)置近場域大小為=3Δ。

圖 12 帶預(yù)制裂紋的受內(nèi)壓圓柱殼模型,mm

根據(jù)模型的對稱性，采用1/2幾何模型，軸向設(shè)置為67個控制點，周向設(shè)置35個控制點，得到2階張量積B樣條曲面。圓柱模型的裂紋擴展路徑見圖13。裂紋沿著與預(yù)置裂紋方向一致的方向擴展，并自始至終沿著與軸線平行的方向向兩側(cè)傳播。在一般殼結(jié)構(gòu)中，所提出的耦合算法在精確幾何表示和斷裂仿真方面的性能均較好。

圖 13 圓柱模型的裂紋擴展路徑

5 結(jié)束語

針對殼結(jié)構(gòu)中的裂紋擴展問題，建立等幾何分析(IGA)與近場動力學(xué)(PD)耦合的分析模型IGA-PD。在計算過程中，將模型域劃分為IGA域與PD域。IGA域采用基于Reissner-Mindlin(R-M)理論的退化殼模型，考慮橫向剪切的影響，保證幾何精確性，避免傳統(tǒng)有限元離散引入的幾何誤差，在法向向量計算等方面更加準(zhǔn)確。PD域采用基于微梁鍵的非連續(xù)伽遼金(DG)有限元殼模型，實現(xiàn)對非連續(xù)場的計算分析。耦合過程基于力/力矩平衡原理，將等幾何剛度矩陣與近場動力學(xué)剛度矩陣耦合。數(shù)值試驗結(jié)果表明，該耦合模型可高效模擬板殼結(jié)構(gòu)中的裂紋擴展，準(zhǔn)確計算裂紋的開裂和交叉等現(xiàn)象。與PD模型相比，本文方法可有效提高計算效率，對于裂紋擴展區(qū)域相對集中的問題，可將計算時間縮短至原耗時的25%。