基于Peridynamics理論的材料可靠性分析數(shù)值方法的優(yōu)化

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

(貴州大學(xué) 機(jī)械工程學(xué)院，貴陽　550025)

?

引用格式：馬騰飛，錢松榮，石宏順，等.基于Peridynamics理論的材料可靠性分析數(shù)值方法的優(yōu)化[J].重慶理工大學(xué)學(xué)報(自然科學(xué)版)，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理論的材料可靠性分析數(shù)值方法的優(yōu)化

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

(貴州大學(xué) 機(jī)械工程學(xué)院，貴陽550025)

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

關(guān)鍵詞：近場動力學(xué)；高斯點(diǎn)積分；優(yōu)化；精度

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

(1)

(2)

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

(3)

s是材料的伸長率，

(4)

則

(5)

其中：s0機(jī)械材料的斷裂極限；K是材料的體積模量?？梢钥闯觯?該方法用積分的方式描述機(jī)械材料的運(yùn)動變形，避免了連續(xù)性假設(shè)和空間微分方程在不連續(xù)狀態(tài)下求解的奇異性，在解決固體力學(xué)斷裂變形中的不連續(xù)問題時具有一定的優(yōu)勢。由于積分在近場動力學(xué)的數(shù)值實(shí)現(xiàn)中有很重要的作用，因此找到合適的積分方法，使得該理論在計算變形破壞有一個合適的誤差精度和一個可接受的計算效率顯得尤為重要。

圖1　xi點(diǎn)的近場區(qū)域

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

1.1離散積分法

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

圖2　CCI方法的晶格的計數(shù)

1.2自適應(yīng)積分法

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

圖3　各節(jié)點(diǎn)積分區(qū)域在k節(jié)點(diǎn)上不同的邊界

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

2固定高斯點(diǎn)方法

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

(6)

表1　位置系數(shù)權(quán)重取值

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

圖4　j節(jié)點(diǎn)的高斯子節(jié)點(diǎn)劃分

在三維坐標(biāo)中子晶格的邊長為

(7)

則子晶格體積為

(8)

將高斯點(diǎn)子節(jié)點(diǎn)應(yīng)用到式(1)得

(9)

(10)

于是，式(9)就變?yōu)?/p>

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

(11)

可以看出：本構(gòu)力的積分公式就轉(zhuǎn)換成了求與一系列高斯點(diǎn)子節(jié)點(diǎn)作用力的和。

3面力密度

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

(12)

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

(13)

4計算實(shí)例

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

由彈塑性力學(xué)應(yīng)力應(yīng)變轉(zhuǎn)換公式[11]得

(14)

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

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

由彈塑性力學(xué)知識知：

(15)

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

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

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

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

5結(jié)論

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

2) 通過改變m值可以看出：對應(yīng)力的計算精度隨Δx的減小而增大。當(dāng)m=8時，計算精度開始趨向穩(wěn)定；但是隨著m值增大，影響域R內(nèi)的相關(guān)節(jié)點(diǎn)增多，增加了計算量，降低了計算效率。

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

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

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

參考文獻(xiàn)：

[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]張民選,羅賢兵.數(shù)值分析[M].南京：南京大學(xué)出版社,2013.

[9]蔡大用.數(shù)值分析與實(shí)驗學(xué)習(xí)指導(dǎo)[M].北京：清華大學(xué)出版社，2002.

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

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

(責(zé)任編輯劉舸)

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

文獻(xiàn)標(biāo)識碼：A

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

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

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

收稿日期：2015-10-23