基于Peridynamics理論的材料可靠性分析數值方法的優化

馬騰飛，錢松榮，石宏順，原群盛

(貴州大學 機械工程學院，貴陽　550025)

?

引用格式：馬騰飛，錢松榮，石宏順，等.基于Peridynamics理論的材料可靠性分析數值方法的優化[J].重慶理工大學學報(自然科學版)，2016(1):32-36.

Citation format：MA Teng-fei，QIAN Song-rong，SHI Hong-shun，et al.Optimization Method for Reliability Analysis of Concrete for the Theory of Peridynamics[J].Journal of Chongqing University of Technology(Natural Science)，2016(1):32-36.

基于Peridynamics理論的材料可靠性分析數值方法的優化

馬騰飛，錢松榮，石宏順，原群盛

(貴州大學 機械工程學院，貴陽550025)

摘要：Peridynamics理論(近場動力學，PD理論)是研究材料破壞變形領域的一門新理論。該理論通過建立相應的數學模型，用積分的方式描述材料的斷裂變形。對PD理論及其傳統數值方法作了簡要的介紹，基于前人針對PD理論提出的積分方法，引入了另一種全新的固定高斯點積分方法。通過對該方法的應用，PD理論數值計算中的效率和精度都有一定改善。

關鍵詞：近場動力學；高斯點積分；優化；精度

Peridynamics理論(PD理論)是2000年美國Sandia國家實驗室的Silling提出的[1-3]。物體中任意兩個物質點之間由于萬有引力的作用能夠產生相互作用力。如圖1所示，該理論基于以上思想，假設有一物質點xi，在xi半徑為δ的影響域R內，物質點xj對xi有相互作用力f，這個力稱為xi點所受到的PD力。通過建立相應的數學模型，把i點在影響域R內所受到的PD力積分起來，再根據牛頓第二定律，給出PD理論的基本運動方程：

(1)

(2)

其中：c是材料的彈性模量，

(3)

s是材料的伸長率，

(4)

則

(5)

其中：s0機械材料的斷裂極限；K是材料的體積模量。可以看出： 該方法用積分的方式描述機械材料的運動變形，避免了連續性假設和空間微分方程在不連續狀態下求解的奇異性，在解決固體力學斷裂變形中的不連續問題時具有一定的優勢。由于積分在近場動力學的數值實現中有很重要的作用，因此找到合適的積分方法，使得該理論在計算變形破壞有一個合適的誤差精度和一個可接受的計算效率顯得尤為重要。

圖1　xi點的近場區域

1PD理論常用的積分方法簡介

1.1離散積分法

由于PD理論積分方法在數值實現中的重要性，Silling提出了把材料離散的積分方法(cubic-cellintergration,CCI)[3-5]。該方法把材料離散成邊長為Δx等大的立方體晶格，如圖2所示。影響域R內對i節點PD力的積分，就轉化為R內各節點對i節點PD力的求和。此方法由于邊界點的計算問題而存在很大誤差。如圖2所示，i節點的實際影響域為以i為圓心以δ半徑為圓形區域，但此方法在計算時，只計算節點在R半徑內的立方體晶格，用實心點表示。如圖2中的晶格2，陰影區域外的一部分也被計算在內，而晶格1本該計算在內的陰影區域內的一部分卻被遺漏掉。

圖2　CCI方法的晶格的計數

1.2自適應積分法

針對Silling CCI方法計算精度低的問題，2010年，B.Kilic提出了一種可以誤差控制的自適應(AI)積分方法，修正了相關點的計算[6-7]，但要對三維坐標系中每個邊界相關點做幾何組態劃分，計算每一個節點的相關體積，導致實現過程很復雜。在該方法中，k既是i的節點又是j的節點，i和j節點影響域與在k節點上的交線分別是AB和CD，那么i和j點有不同的積分邊界，在AI方法中，此處k就有2組不同的梯形點，如圖3所示。另外，PD力積分需要積分點的最近位移信息，然而由運動方程決定的大量梯形積分點的位移需要每一步單獨計算一次，這都使得該方法計算效率較低。

圖3　各節點積分區域在k節點上不同的邊界

為了避免不同積分極限帶來的麻煩，且在可以接受的效率內保證精度，本研究引入固定高斯點計算方法。

2固定高斯點方法

高斯(Gauss)積分法是數值分析積分中常用的一種近似求積方法[8-9]。Gauss積分公式為

(6)

表1　位置系數權重取值

正確找到高斯點位置后，在坐標中進一步劃分相關節點，使得一個高斯點都位于一個子節點的中心。如圖4所示，“×”代表高斯點，圖4(a)中將一個晶格劃分為4個子晶格，圖4(b)中將一個晶格劃分為9個子晶格。

圖4　j節點的高斯子節點劃分

在三維坐標中子晶格的邊長為

(7)

則子晶格體積為

(8)

將高斯點子節點應用到式(1)得

(9)

(10)

于是，式(9)就變為

f(ug-ui,xg-xi)Vg+bi

(11)

可以看出：本構力的積分公式就轉換成了求與一系列高斯點子節點作用力的和。

3面力密度

假設有一均勻變形的物體R，在R內取一點i，讓通過i的單位矢量n的法平面將R分成兩部分R+和R-。令：

(12)

則i的n方向上的點的面力密度τ(i,n)為

(13)

4計算實例

現假設有一個均勻的各項同性的微彈性脆性材料的彈性模量為1.0×105N/mm2，有效積分區域為δ，泊松比為0.25，極限拉伸系數[10]s0=0.001 2。這個微彈材料受到一個Y方向上的應變矩陣：

由彈塑性力學應力應變轉換公式[11]得

(14)

此處取E=1.0×105N/mm2。由彈塑性力學知識可得其應力矩陣如下：

根據式(11)將此微彈材料進行數學建模，以i點為所求面力的中心點，i的影響域為δ，將δ均分m段，令m=12，即Δx=δ/12。取式(6)中n=2，即在三維情況下把1個晶格劃分為8個高斯點子晶格，那么R+區域就離散為很多含有一個高斯點xg的小晶格，相對于i點的體力積分就轉換為許多xg晶格相對于中心點i的體力的求和。

由彈塑性力學知識知：

(15)

其中u1,u2,u3分別為點xg相對于i點的3個方向上的位移，即式(11)中的ug～ui，x,y,z分別為3個方向上xg點相對i點的相對位置，為式(11)中的xg～xi。

依次改變Δx的大小，編程計算得到應力張量，與理論值對比，得出誤差曲線隨m值的變化，如圖4所示。

圖5　精度與Δx值之間關系(n=2)

取式(6)中n=3，即在三維坐標中把1個晶格劃分為27個高斯子晶格，仍取Δx=δ/12，編程計算得到相應的應力張量為

5結論

1) 將通過編程計算得到的應力矩陣與實際的應力矩陣對比可知：PD理論的固定高斯點方法計算所得的結果與實際數據相符(誤差為4.09%)，證明該方法是可行的。

2) 通過改變m值可以看出：對應力的計算精度隨Δx的減小而增大。當m=8時，計算精度開始趨向穩定；但是隨著m值增大，影響域R內的相關節點增多，增加了計算量，降低了計算效率。

3) 當高斯點精細時，隨著子節點的增多，計算精度得到明顯提高(當n=3，m=12時，誤差為1.05%)。

4) 當m值相同時，固定高斯點方法的計算精度高于Silling提出的CCI方法的精度[3]，但計算時間略長。對比AI方法的實現過程及精度[6]可知：固定高斯點方法的精度略低于AI方法,但實現過程大為簡化，計算效率有明顯提高。

通過以上的計算可以看出：固定高斯點積分法是PD理論運動方程數值求解方法的改進，為PD理論的數值積分的實現，為對精度和效率有不同需求的情況提供了一條切實可行的途徑。

參考文獻：

[1]SILLING S A. Reformulation of elasticity theory for discontinuities and long-range forces[J].Journal of the Mechanics and Physics of Solids,2000,48(1)：175-209.

[2]SILLING S A,ASKARI E.Peridynamic modeling of impact damage[C]//Proceeding of the ASME/JSME Pressure Vessels and Piping Conference.San Diego:[s.n.],2000:197-205.

[3]SILLING S A,ASKARI E.A meshfree method based on the peridynamic model of solid mechanics[J].Computers and Structures,2005,83:1526-1535.

[4]SILLING S A. Dynamic fracture modeling with a meshfree peridynamic code.Second MIT Conference on Computional Fluid and Solid Mechanics.Massachusettes：MIT,2003:641-644.

[5]MACEK R W,SILLING S A.Peridynamics via finite element analysis[J].Finite Elements in Analysis and Design,2007,43(15):1169-1178.

[6]YU K,XIN X J,LEASE K B.A new method of adaptive integration with error control for bond-based peridynamics[C]//Proceeding of the World Congress on Engineering and Computer Science.San Francisco:[s.n.],2010:1041-1046.

[7]YU K,XIN X J，LEASE K B.A new adaptive integration method for the peridynamic theory[J].Modeling and Simulation in Materials Science and Engineering,2011,19:45003.

[8]張民選,羅賢兵.數值分析[M].南京：南京大學出版社,2013.

[9]蔡大用.數值分析與實驗學習指導[M].北京：清華大學出版社，2002.

[10]沈峰，章青，黃丹，等.基于近場動力學理論的混凝土軸拉破壞過程模擬[J].計算力學學報,.2013,30(6):79-83.

[11]陳明祥.彈塑性力學[M].武漢：科學出版社，2010.

(責任編輯劉舸)

Optimization Method for Reliability Analysis of

Concrete for the Theory of Peridynamics

MA Teng-fei，QIAN Song-rong，SHI Hong-shun，YUAN Qun-sheng

(School of Mechanical Engineering, Guizhou University, Guiyang 550025，China)

Abstract:Peridynamics (PD) is an emerging new theory used in the research of material damage. This theory describes the fracture and deformation of material by establishing corresponding mathematical model and using integral method. This paper simply introduced the theory of peridynamics and traditional numerical method and a new integration method with fixed Gaussian points was introduced based on previous methods. The efficiency and accuracy of the numerical calculation of PD theory also have relative improvement through the application of the new method.

Key words:peridynamics; Gauss integral; optimization; accuracy

文章編號：1674-8425(2016)01-0032-05

中圖分類號：TH115

文獻標識碼：A

doi:10.3969/j.issn.1674-8425(z).2016.01.006

作者簡介：馬騰飛(1989—)，男，河南商丘人，碩士研究生，主要從事機械材料的可靠性分析研究。

基金項目：貴州省國際合作項目(黔科合外G字[2013]7006號)；貴州省聯合 (黔科合LH字[2014]7624號)；貴州省留學人員科技活動項目 (黔人項目資助合同[2014]13號)

收稿日期：2015-10-23