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

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

?

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

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

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

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

(1)

(2)

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

(3)

s是材料的伸長(zhǎng)率，

(4)

則

(5)

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

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

1PD理論常用的積分方法簡(jiǎn)介

1.1離散積分法

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

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

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

針對(duì)Silling CCI方法計(jì)算精度低的問題，2010年，B.Kilic提出了一種可以誤差控制的自適應(yīng)(AI)積分方法，修正了相關(guān)點(diǎn)的計(jì)算[6-7]，但要對(duì)三維坐標(biāo)系中每個(gè)邊界相關(guān)點(diǎn)做幾何組態(tài)劃分，計(jì)算每一個(gè)節(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)動(dòng)方程決定的大量梯形積分點(diǎn)的位移需要每一步單獨(dú)計(jì)算一次，這都使得該方法計(jì)算效率較低。

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

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

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

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

(6)

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

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

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

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

(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計(jì)算實(shí)例

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

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

(14)

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

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

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

(15)

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

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

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

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

5結(jié)論

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

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

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

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

通過以上的計(jì)算可以看出：固定高斯點(diǎn)積分法是PD理論運(yùn)動(dòng)方程數(shù)值求解方法的改進(jìn)，為PD理論的數(shù)值積分的實(shí)現(xiàn)，為對(duì)精度和效率有不同需求的情況提供了一條切實(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í)驗(yàn)學(xué)習(xí)指導(dǎo)[M].北京：清華大學(xué)出版社，2002.

[10]沈峰，章青，黃丹，等.基于近場(chǎng)動(dòng)力學(xué)理論的混凝土軸拉破壞過程模擬[J].計(jì)算力學(xué)學(xué)報(bào),.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

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

中圖分類號(hào)：TH115

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

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

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

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

收稿日期：2015-10-23