基于Kriging插值無(wú)網(wǎng)格法的動(dòng)力彈塑性分析

熊勇剛 ，田萬(wàn)鵬，陳科良，陳云飛 ，杜民獻(xiàn)，吳吉平

(1. 湖南工業(yè)大學(xué) 機(jī)械工程學(xué)院，湖南 株洲，412007；2. 東南大學(xué) 機(jī)械工程學(xué)院，江蘇 南京，210000)

無(wú)網(wǎng)格方法近年來(lái)得到了較大發(fā)展，是目前科學(xué)和工程計(jì)算方法研究的熱點(diǎn)之一[1?3]。與有限元法不同，無(wú)網(wǎng)格方法的近似函數(shù)是建立在一系列離散點(diǎn)上的，可以部分或徹底消除網(wǎng)格，從而有效地避免了由于網(wǎng)格的依賴(lài)性而帶來(lái)的種種困難。目前，發(fā)展的無(wú)網(wǎng)格方法主要有無(wú)單元 Galerkin法[4?5]、無(wú)網(wǎng)格局部Petrov-Galerkin法[6]、光滑粒子流體動(dòng)力學(xué)法[7?8]、再生核粒子法[9]、復(fù)變量無(wú)網(wǎng)格方法[10]、邊界無(wú)單元法[11?12]、無(wú)網(wǎng)格流形方法[13]以及自然單元法[14]等，其中，無(wú)單元Galerkin法是發(fā)展良好和最具有應(yīng)用價(jià)值的無(wú)網(wǎng)格方法之一[1]，但該法的最大缺點(diǎn)是采用移動(dòng)最小二乘法構(gòu)造的形函數(shù)不滿(mǎn)足Kronecker delta函數(shù)性質(zhì)，因此，本質(zhì)邊界條件的施加非常困難，需要采用特別技術(shù)進(jìn)行處理。Kriging插值無(wú)網(wǎng)格法是在構(gòu)造形函數(shù)時(shí)采用滑動(dòng) Kriging插值代替移動(dòng)最小二乘法而發(fā)展起來(lái)的一種改進(jìn)的無(wú)單元 Galerkin法。滑動(dòng)Kriging插值法構(gòu)造近似函數(shù)不僅數(shù)值實(shí)現(xiàn)簡(jiǎn)單、計(jì)算量小，而且滿(mǎn)足Kronecker delta函數(shù)性質(zhì)，從而可以非常方便地施加本質(zhì)邊界條件[15?16]。為此，本文作者將 Kriging插值無(wú)網(wǎng)格法應(yīng)用于動(dòng)力彈塑性問(wèn)題的求解計(jì)算。首先詳細(xì)推導(dǎo) Kriging插值無(wú)網(wǎng)格法在動(dòng)力彈塑性問(wèn)題中的理論公式，然后，給出基于 Kriging插值無(wú)網(wǎng)格法的動(dòng)力彈塑性分析的數(shù)值實(shí)現(xiàn)方法和具體算例。

1 滑動(dòng)Kriging插值

與移動(dòng)最小二乘法相似，假設(shè)任意計(jì)算點(diǎn)x的影響域內(nèi)包含n個(gè)節(jié)點(diǎn)。滑動(dòng)Kriging插值的逼近函數(shù)采用如下線(xiàn)性回歸模型與偏差之和的模式：

式中：pT(x)為多項(xiàng)式基函數(shù)；a為待求常系數(shù)向量；z(x)為一個(gè)均值為0、方差為σ2、協(xié)方差不為0時(shí)的局部離差。對(duì)于二維情況，p(x)可選取如下二次基：

用式(1)對(duì)n個(gè)樣本點(diǎn)進(jìn)行插值，則z(x)的協(xié)方差可表示為

θ＞0，為模型參數(shù)；為給定2點(diǎn)xi和xj之間的相關(guān)距離。

計(jì)算點(diǎn)x與所給定的n個(gè)節(jié)點(diǎn)之間的相關(guān)函數(shù)向量為

在給定的n個(gè)節(jié)點(diǎn)處的函數(shù)值已知時(shí)，

式中：u為n個(gè)節(jié)點(diǎn)處場(chǎng)函數(shù)的集合。當(dāng)采用線(xiàn)性回歸模型式(1)表達(dá)u時(shí)，有

式中：Pm與Z分別為給定點(diǎn)上基函數(shù)值所形成的矩陣及線(xiàn)性逼近的誤差向量，且有

對(duì)計(jì)算點(diǎn)x，采用已知節(jié)點(diǎn)的函數(shù)值對(duì)場(chǎng)函數(shù)u(x)進(jìn)行線(xiàn)性逼近，有

于是，式(1)與式(10)之間的誤差函數(shù)為

要保證無(wú)偏估計(jì)，必須滿(mǎn)足如下約束條件：

式(11)中誤差函數(shù)的均方差為

為了確定最優(yōu)的線(xiàn)性估計(jì)，將式(13)作為目標(biāo)函數(shù)，式(12)作為約束條件，進(jìn)而根據(jù)Lagrange不定乘子法可以將其轉(zhuǎn)化為無(wú)條件的極小值問(wèn)題。通過(guò)對(duì)向量λT(x)求偏導(dǎo)數(shù)為0得到一組方程，求解此方程組可得到最優(yōu)線(xiàn)性無(wú)偏估計(jì)條件下的插值函數(shù)[15?16]：

式中：

同時(shí)，為了便于表述，定義

式中：I為n×n階單位矩陣。于是，式(14)可改寫(xiě)為

AjI為矩陣A的第j行、第I列元素；BkI為矩陣B的第k行、第I列元素。

形函數(shù)φI(x)關(guān)于x和y的導(dǎo)函數(shù)分別為

2 動(dòng)力彈塑性分析的Kriging插值法

2.1 控制方程

考慮二維動(dòng)力彈塑性問(wèn)題，其計(jì)算域?yàn)?，邊界為控制方程為：

式中：ρ為質(zhì)量密度；為節(jié)點(diǎn)加速度，且有=?2u/?t2；ui為節(jié)點(diǎn)位移；t為時(shí)間；為應(yīng)力；為應(yīng)變；bi為體力；Γu和Γt分別為位移和面力已知的邊界。采用von Mises屈服準(zhǔn)則、關(guān)聯(lián)流動(dòng)準(zhǔn)則和各向同性硬化準(zhǔn)則時(shí)，式(22)中的彈塑性

其中：

σeq和H分別為等效應(yīng)力和塑性模量；G和v分別為剪切模量和泊松比。

應(yīng)力和位移滿(mǎn)足如下邊界條件

和初始條件

式中：nj(x)為Γ上點(diǎn)x處的外法線(xiàn)方向余弦；和分別為已知的位移和面力分量；u0和v0分別為初始位移和初始速度。

2.2 空間離散

平衡方程式(21)的弱形式可以通過(guò)加權(quán)殘值法獲得，即

式中：wi為權(quán)函數(shù)。對(duì)式(29)進(jìn)行分部積分，并考慮到自然邊界條件式(27)，有

由于只對(duì)空間域進(jìn)行離散，求解域?內(nèi)的試函數(shù)可以由式(18)表示為

將空間離散后的位移表達(dá)式(31)代入式(30)，且令權(quán)函數(shù)等于形函數(shù)，最終可得彈性動(dòng)力學(xué)的離散方程為

式中：M，K和f(t)分別為質(zhì)量矩陣、剛度矩陣與載荷向量；

對(duì)于動(dòng)力彈塑性問(wèn)題，式(32)可以改寫(xiě)為

節(jié)點(diǎn)內(nèi)力向量P(σ(t))為

其中：σ(t)為應(yīng)力向量。

2.3 求解方案

本文采用傳統(tǒng)的 Newmark方法對(duì)時(shí)間域進(jìn)行離散。Newmark法于1959年由線(xiàn)性加速度法延伸推導(dǎo)而得，是一種應(yīng)用相當(dāng)廣泛的直接積分法。由于塑性變形的非線(xiàn)性特性，在每個(gè)時(shí)間步都必須進(jìn)行迭代計(jì)算。在t和t+Δt時(shí)間步內(nèi)的迭代求解過(guò)程如下[17?19]：

(1) 令迭代計(jì)算變量k=0；

(2) 在開(kāi)始預(yù)估階段，令

(3) 利用下述方程計(jì)算殘余力：

(4) 如有需要，應(yīng)用下式形成等效剛度矩陣：

否則，就應(yīng)用前面已計(jì)算的K*。

(5) 求解方程：

(6) 進(jìn)入修正階段，令

式中：α和δ是按積分精度和穩(wěn)定性要求決定的參數(shù)。當(dāng)δ≥0.5和α≥0.25(0.5+δ)2時(shí)算法是無(wú)條件穩(wěn)定的，即時(shí)間步長(zhǎng)Δt不影響解的穩(wěn)定性。

(7) 迭代收斂檢查。對(duì)于給定的收斂誤差εR，若||Δu(k)||≤εR，則結(jié)束本時(shí)間步內(nèi)的迭代循環(huán)，否則，令k=k+1并轉(zhuǎn)到第3步。

(8) 為了使下一個(gè)時(shí)間步使用，令

給定初始位移u0和初始速度v0后，可以從下式求出初始加速度：

3 數(shù)值算例

為了驗(yàn)證以上提出的數(shù)值方法的有效性，本文分析2個(gè)算例。在滑動(dòng)Kriging插值中，本文選用二次基函數(shù)，并且計(jì)算點(diǎn)x的影響域取為以點(diǎn)x為圓心的圓形域，半徑dmax=scale×s[6](其中s[6]為計(jì)算點(diǎn)x到距其最近的第 6個(gè)節(jié)點(diǎn)之間的距離)，系數(shù)scale取為2.5。

算例1 受內(nèi)壓作用的厚壁圓筒。

1個(gè)內(nèi)徑a=1 m，外徑b=2 m的厚壁圓筒，受到p=1.85×104Pa的突加內(nèi)壓作用，如圖1所示。材料為理想彈塑性，并服從von Mises屈服準(zhǔn)則。材料的彈性模量E=2.1×107Pa，泊松比ν=0.3，屈服應(yīng)力σy=3.55×104Pa，質(zhì)量密度ρ=7.85 kg/m3。該問(wèn)題為平面應(yīng)變狀態(tài)下的1個(gè)經(jīng)典算例。

圖1 受內(nèi)壓作用的厚壁圓筒Fig. 1 A thick-walled cylinder subjected to internal pressure

由于對(duì)稱(chēng)性，可取結(jié)構(gòu)的1/4作為計(jì)算模型。計(jì)算分析中采用圖 1(a)所示的節(jié)點(diǎn)布置方案，且時(shí)間步長(zhǎng)取為Δt=1.0×10?5s。圖2所示為內(nèi)表面的徑向位移隨時(shí)間的變化曲線(xiàn)。為了進(jìn)行對(duì)比，圖2還給出了有限元軟件ABAQUS的計(jì)算結(jié)果。從圖2可以看出：本文的計(jì)算結(jié)果與 ABAQUS的計(jì)算結(jié)果很吻合，說(shuō)明本文推導(dǎo)的動(dòng)力彈塑性分析的 Kriging插值無(wú)網(wǎng)格法及其離散形式是正確的，方法是可行的。

圖2 內(nèi)表面的徑向位移Fig. 2 Radial displacement history at inner surface

算例2 受突加載荷的懸臂梁。

如圖3所示為長(zhǎng)L=1.2 m、高H=0.4 m的懸臂梁，在平面應(yīng)力狀態(tài)下端部受到突加剪切荷載p=2.5 kPa的作用。材料為理想彈塑性，并服從von Mises屈服準(zhǔn)則。材料彈性模量E=2.1×107Pa，泊松比ν=0.3，屈服應(yīng)力σy=4.0×104Pa，質(zhì)量密度ρ=7.85 kg/m3。

在計(jì)算時(shí)，把區(qū)域均勻地劃分為25×7個(gè)節(jié)點(diǎn)的網(wǎng)格，且時(shí)間步長(zhǎng)取為Δt=1.0×10?4s。為了檢驗(yàn)方法的有效性，本文使用有限元軟件 ABAQUS的計(jì)算結(jié)果進(jìn)行對(duì)比。圖4所示為懸臂梁右端中點(diǎn)A的豎向位移隨時(shí)間的變化曲線(xiàn)。從圖4可以看出：本文的計(jì)算結(jié)果與 ABAQUS的計(jì)算結(jié)果很吻合，從而進(jìn)一步驗(yàn)證了本文方法的有效性。

圖3 受突加載荷的懸臂梁Fig. 3 Cantilever beam under Heaviside loading

圖4 懸臂梁中A點(diǎn)的豎向位移Fig. 4 Vertical displacement of point A in beam

4 結(jié)論

(1) 采用Kriging插值無(wú)網(wǎng)格法，建立了進(jìn)行結(jié)構(gòu)彈塑性動(dòng)力響應(yīng)分析的一種新方法，推導(dǎo)了相應(yīng)的計(jì)算公式，提出了具體的數(shù)值實(shí)現(xiàn)方法。數(shù)值計(jì)算結(jié)果表明，采用 Kriging插值無(wú)網(wǎng)格法進(jìn)行動(dòng)力彈塑性分析是可行的，它具有穩(wěn)定性好和精度高等優(yōu)點(diǎn)。

(2) 采用Kriging插值無(wú)網(wǎng)格法進(jìn)行動(dòng)力彈塑性分析只需要離散節(jié)點(diǎn)的信息，從而減少了前處理的工作量。尤為重要的是，滑動(dòng) Kriging插值法的形函數(shù)滿(mǎn)足Kronecker delta函數(shù)性質(zhì)，因而能夠直接、準(zhǔn)確地施加本質(zhì)邊界條件。此外，形函數(shù)導(dǎo)數(shù)的計(jì)算不涉及逆矩陣的求導(dǎo)，因而計(jì)算效率較高。

[1]Belytschko T, Krongauz Y, Organ D, et al. Meshless methods:An overview and recent developments[J]. Computer Methods in Applied Mechanics and Engineering, 1996, 139: 3?47.

[2]Li S F, Liu W K. Meshfree and particle methods and their applications[J]. Applied Mechanics Review, 2002, 55: 1?34.

[3]宋康祖, 陸明萬(wàn), 張雄. 固體力學(xué)中的無(wú)網(wǎng)格方法[J]. 力學(xué)進(jìn)展, 2000, 30(3): 55?65.SONG Kangzu, LU Mingwan, ZHANG Xiong. Meshless method for solid mechanics[J]. Advance in Mechanics, 2000,30(3): 55?65.

[4]Belytschko T, Lu Y Y, Gu L. Element-free Galerkin method[J].International Journal for Numerical Methods in Engineering,1994, 37: 229?256.

[5]王東東, 李凌, 張燦輝. 穩(wěn)定節(jié)點(diǎn)積分伽遼金無(wú)網(wǎng)格法的應(yīng)力計(jì)算方法研究[J]. 固體力學(xué)學(xué)報(bào), 2009, 30(6): 586?591.WANG Dongdong, LI Ling, ZHANG Canhui. On stress evaluation in Galerkin meshfree methods with stabilized conforming nodal integration[J]. Chinese Journal of Solid Mechanics, 2009, 30(6): 586?591.

[6]Atluri S N, Zhu T L. A new meshless local Petrov-Galerkin(MLPG) approach in computational mechanics[J]. Computational Mechanics, 1998, 22(2): 117?127.

[7]Gingold R A, Monaghan J J. Smoothed particle hydrodynamics:theory and application to non-spherical stars[J]. Monthly Notices of the Royal Astronomical Society, 1977, 181: 375?389.

[8]繆吉倫, 陳景秋, 張永祥. 基于 SPH 方法的立面二維涌浪數(shù)值模擬[J]. 中南大學(xué)學(xué)報(bào)(自然科學(xué)版), 2012, 43(8):3244?3249.MIAO Jilun, CHEN Jingqiu, ZHANG Yongxiang.Two-dimension numerical simulation of impulsive waves by SPH method[J]. Journal of Central South University (Science and Technology), 2012, 43(8): 3244?3249.

[9]Liu W K, Jun S, Zhang Y F. Reproducing kernel particle methods[J]. International Journal for Numerical Methods in Fluids, 1995, 20: 1081?1106.

[10]程玉民, 李九紅. 彈性力學(xué)的復(fù)變量無(wú)網(wǎng)格方法[J]. 物理學(xué)報(bào), 2005, 54(10): 4463?4471.CHENG Yumin, LI Jiuhong. A meshless method with complex variables for elasticity[J]. Acta Physica Sinica, 2005, 54(10):4463?4471.

[11]程玉民, 陳美娟. 彈性力學(xué)的一種邊界無(wú)單元法[J]. 力學(xué)學(xué)報(bào), 2003, 35(2): 181?186.CHENG Yumin, CHEN Meijuan. A boundary element-free method for linear elasticity[J]. Chinese Journal of Theoretical and Applied Mechanics, 2003, 35(2): 181?186.

[12]程玉民, 彭妙娟. 彈性動(dòng)力學(xué)的邊界無(wú)單元法[J]. 中國(guó)科學(xué),2005, 35(4): 435?448.CHENG Yumin, PENG Miaojuan. Boundary element-free method for elastodynamics[J]. Science in China, 2005, 35(4):435?448.

[13]李樹(shù)忱, 李術(shù)才, 朱維申. 彈性動(dòng)力學(xué)的無(wú)網(wǎng)格流形方法[J].巖石力學(xué)與工程學(xué)報(bào), 2006, 25(1): 141?145.LI Shuchen, LI Shucai, ZHU Weishen. Study on meshless manifold method in elastodynamics[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(1): 141?145.

[14]Sukumar N, Moran B, Belytschko T. The natural elements method in solid mechanics[J]. International Journal for Numerical Methods in Engineering, 1998, 43: 839?887.

[15]Gu L. Moving kriging interpolation and element-free Galerkin method[J]. International Journal for Numerical Methods in Engineering, 2003, 56: 1?11.

[16]Dai K Y, Liu G R, Lim K M, et al. Comparison between the radial point interpolation and the kriging interpolation used in meshfree methods[J]. Computational Mechanics, 2003, 32:60?70.

[17]Soares D Jr, Sladek J, Sladek V. Nonlinear dynamic analyses by meshless local Petrov-Galerkin formulations[J]. International Journal for Numerical Methods in Engineering, 2010, 81:1687?1699.

[18]陳莘莘, 李慶華, 劉應(yīng)華, 等. 動(dòng)力彈塑性分析的無(wú)網(wǎng)格自然單元法[J]. 固體力學(xué)學(xué)報(bào), 2011, 32(5): 493?499.CHEN Shenshen, LI Qinghua, LIU Yinghua, et al. Application of meshless natural element method to dynamic elastoplastic analysis[J]. Chinese Journal of Solid Mechanics, 2011, 32(5):493?499.

[19]Owen D J R, Hinton E. Finite element in plasticity: Theory and practice[M]. Swansea: Pineridge Press, 1980: 431?463.