房霆宸， 李 武

(1.上海大學(xué) 土木工程系， 上海 200072；2.中交第三航務(wù)工程勘察設(shè)計(jì)院有限公司， 上海 200032)

無(wú)網(wǎng)格法在動(dòng)力學(xué)中的理論研究和應(yīng)用已吸引了許多國(guó)內(nèi)外學(xué)者，是計(jì)算力學(xué)中的一個(gè)研究熱點(diǎn)[1]。自然單元法[2～6](Natural Element Method，NEM)是最近出現(xiàn)的一種無(wú)網(wǎng)格方法，自1995年Braun和Sambridge 在Nature上發(fā)表的《A numerical method for solving partial differ-ential equations on highly irregular evolving grids》一文提出自然單元法以來(lái)，以Voronoi圖為幾何基礎(chǔ)的數(shù)值方法在國(guó)內(nèi)外得到極大的關(guān)注.許多學(xué)者開(kāi)展研究，并在不同的工程技術(shù)領(lǐng)域得到一定的應(yīng)用。1998年，Sukumar N[3]將自然單元法成功地應(yīng)用于求解固體力學(xué)中的橢圓型偏微分方程，并且構(gòu)造了C1自然相鄰節(jié)點(diǎn)插值成功求解了橢圓型四階偏微分方程。在此階段，NEM形函數(shù)的構(gòu)造采用的是Sibson插值方法。2001年，Sukumar N[5]采用了Non-Sibsonian插值方法構(gòu)造NEM形函數(shù)，從而使得NEM方法更為方便實(shí)用。Cueto E[6]基于α-shape的概念實(shí)現(xiàn)了三維分析。蔡永昌、朱合華等[7]研究了采用局部Petrov-Galerkin方法獲得整體系統(tǒng)平衡的控制方程。盧波、葛修潤(rùn)等[8，9]研究了弱形式的數(shù)值積分方案，對(duì)自然單元法的原始應(yīng)力數(shù)值解的精度進(jìn)行了改善提高。在應(yīng)用方面，Martniez M A[2，10]采用自然單元法求解了流體力學(xué)問(wèn)題，朱合華、楊寶紅[11]采用自然單元法求解了彈塑性問(wèn)題。

無(wú)網(wǎng)格方法在動(dòng)力學(xué)中得到廣泛應(yīng)用，自然單元法作為無(wú)網(wǎng)格方法的一種也應(yīng)該被應(yīng)用到動(dòng)力學(xué)中。因此，本文詳細(xì)推導(dǎo)和討論自然單元法在動(dòng)力學(xué)問(wèn)題中離散方程和求解方法。

1 位移函數(shù)

1.1 近似位移函數(shù)

仿照二維Voronoi圖的構(gòu)建理論建立空間8結(jié)點(diǎn)的Voronoi圖結(jié)構(gòu)。本文采用立方體的頂點(diǎn)作為結(jié)點(diǎn)對(duì)空間進(jìn)行劃分，得到8結(jié)點(diǎn)的Voronoi結(jié)構(gòu)，如圖1所示。在圖1中放入待插值點(diǎn)X，并對(duì)8個(gè)結(jié)點(diǎn)Voronoi圖進(jìn)行二次劃分，得到待插值點(diǎn)X的二階Voronoi結(jié)構(gòu)的1/8部分，如圖2所示。根據(jù)空球規(guī)則尋找自然相鄰結(jié)點(diǎn)——若四面體的外接球包含待插值點(diǎn)，則該四面體的4個(gè)頂點(diǎn)即為待插值點(diǎn)的自然相鄰結(jié)點(diǎn)。本文中采用的是立方體，再有本文四面體的外接球等價(jià)立方體的外接球，因此8個(gè)頂點(diǎn)既是待插值點(diǎn)X的自然相鄰結(jié)點(diǎn)。

圖1 8結(jié)點(diǎn)一階Voronoi圖

標(biāo)號(hào)為1～8的8個(gè)結(jié)點(diǎn)構(gòu)成了待插值點(diǎn)X的自然相鄰結(jié)點(diǎn)，待插值點(diǎn)與自然相鄰點(diǎn)一起形成Voronoi圖稱(chēng)為待插值點(diǎn)的二次Voronoi結(jié)構(gòu)。因?yàn)榭臻g結(jié)構(gòu)復(fù)雜，所以采用1/8結(jié)構(gòu)進(jìn)行說(shuō)明推導(dǎo)過(guò)程。

圖2 8結(jié)點(diǎn)二階Voronoi圖

在圖1、圖2中1～8為空間結(jié)點(diǎn)，A～H為7結(jié)點(diǎn)的Voronoi圖頂點(diǎn)；在圖1中X為插值點(diǎn)，與A重合，X與8個(gè)結(jié)點(diǎn)形成二階Voronoi圖結(jié)構(gòu)中的1/8部分，其中abc為X點(diǎn)二階Voronoi圖中3個(gè)頂點(diǎn)，取出與7結(jié)點(diǎn)相關(guān)部分的四面體Aabc，四面體Aabc為X的二階Voronoi圖與7結(jié)點(diǎn)一階Voronoi圖的重疊部分，而且直線A7垂直平面abc。根據(jù)圖形對(duì)稱(chēng)性知道，其它7個(gè)結(jié)點(diǎn)的重疊部分與7結(jié)點(diǎn)相似，X的二階Voronoi既是由8個(gè)四面體組成多面體。本文采用多面體相關(guān)邊長(zhǎng)建立函數(shù)。下面以7結(jié)點(diǎn)為例說(shuō)明構(gòu)建位移插值函數(shù)過(guò)程。

(1)

(2)

式中，S7(x)=ac+ab+bc是與結(jié)點(diǎn)7關(guān)聯(lián)的Voronoi邊的長(zhǎng)度之和，h7(x)=A7是插值點(diǎn)x到結(jié)點(diǎn)7的Voronoi邊組成平面的垂直距離(圖2)。

把式(1)、式(2)中的7擴(kuò)展為域內(nèi)任意鄰接點(diǎn)I的插值函數(shù)：

(3)

(4)

式中，SI(x)是與結(jié)點(diǎn)I關(guān)聯(lián)的Voronoi邊的長(zhǎng)度之和，hI(x)是插值點(diǎn)x到結(jié)點(diǎn)I的Voronoi邊組成平面的垂直距離。

1.2 三維位移插值函數(shù)的性質(zhì)

由方程(3)的形函數(shù)構(gòu)造過(guò)程可以很容易證明：

0≤ΦI(x)≤1

(5)

在圖2中，可以注意到如果點(diǎn)x與任何結(jié)點(diǎn)重合，例如與結(jié)點(diǎn)7重合，則Φ7(x)=1，ΦI(x)=0(I≠7)。因此，該形函數(shù)與有限元法形函數(shù)一樣，滿足Kronecker條件：

ΦI(xJ)=δIJ

(6)

由式(3)可知，與有限元法一樣，自然鄰接形函數(shù)還滿足單位分解條件和線性連續(xù)條件：

(7)

(8)

此外，自然鄰接點(diǎn)的形函數(shù)除了在結(jié)點(diǎn)處C0連續(xù)外，在其它地方C∞連續(xù)。

1.3 三維位移插值函數(shù)的導(dǎo)數(shù)

根據(jù)插值函數(shù)式(3)、式(4)求導(dǎo)相應(yīng)導(dǎo)數(shù)如下：

(9)

(10)

2 控制方程

平衡方程

σij,j(x,t)+fi(x,t)-ρ(x)ui,tt(x,t)-

μui,t(x,t)=0 (在Ω域內(nèi))

(11)

幾何方程

(12)

物理方程

σij(x,t)=Dijklεkl(x,t) (在Ω域內(nèi))

(13)

邊界條件

(14)

(15)

初始條件

(16)

(17)

平衡方程(11)式及力的邊界條件(15)式的等效積分形式的伽遼金表示法如下：

(18)

δui(x,t)ρui,tt(x,t)+δui(x,t)μui,t(x,t))dΩ

(19)

式中，δui(x,t)是任意變化的空間位移增量。

采用自然單元法對(duì)空間域進(jìn)行離散，所以得到離散域內(nèi)任意點(diǎn)x的位移u，v，w插值為：

(20)

(21)

(22)

式中，Φ(x)為結(jié)點(diǎn)插值的形函數(shù)，x為結(jié)點(diǎn)空間幾何坐標(biāo)。

將空間離散后的位移表達(dá)式(20)至式(22)(u(x,t)=u1(x,t),v(x,t)=u2(x,t)，w(x,t)=u3(x,t))代入(19)式得到系統(tǒng)求解方程如下：

(23)

式中：

(24)

(25)

(26)

(27)

式中，M，C，K，F(xiàn)分別為質(zhì)量矩陣、阻尼矩陣、剛度矩陣及荷載向量。

3 算例分析

通過(guò)以下算例，將本文方法與自然單元法(NEM)進(jìn)行了對(duì)比。

為了對(duì)誤差進(jìn)行比較分析，定義位移誤差范數(shù)Lu：

Lu=

(28)

3.1 分片試驗(yàn)

驗(yàn)算基于non-Sibsonian插值的三維自然單元法的收斂性，引入有限元中位移分片試驗(yàn)來(lái)驗(yàn)證本文程序的正確性。被檢驗(yàn)Ω是立方體，如圖3所示，采用均勻布點(diǎn)和非均勻布點(diǎn)兩種離散方法，結(jié)點(diǎn)數(shù)為216、238，如圖4、5所示。在立方體表面結(jié)點(diǎn)賦予該結(jié)點(diǎn)坐標(biāo)作為給定位移，比較該立方體內(nèi)部結(jié)點(diǎn)位移值ui與其坐標(biāo)值。并在被檢驗(yàn)Ω是立方體中任意選取8個(gè)節(jié)點(diǎn)，比較該立方體內(nèi)部結(jié)點(diǎn)位移值ui與其坐標(biāo)值，如表1。

圖3 立方體

圖4 均勻布點(diǎn)

圖5 非均勻布點(diǎn)

編號(hào)a位移值x方向a坐標(biāo)值x方向b位移值x方向b坐標(biāo)值x方向10.20280.20280.101970.1019720.20740.20740.134740.1347430.40490.40490.242360.2423640.40730.40730.496730.4967350.60100.60100.596380.5963860.60270.60100.608130.6081370.80150.60100.727150.7271580.80500.60100.843280.84328

由表1可知本文三維自然單元法的算法可以精確通過(guò)分片試驗(yàn)，因此該算法解的收斂性通過(guò)分片試驗(yàn)的檢驗(yàn)，說(shuō)明本文算法是正確的、可行的。

3.2 懸臂梁端部受集中靜載荷作用

三維彈性懸臂梁純彎曲小變形狀態(tài)的撓度解析解為：

(29)

(30)

(31)

式中，p為集中荷載，E為彈性模量，v為剪切模量，I為矩形截面慣性矩。

懸臂梁尺寸如下：l=4000，b=h=400，如圖6。右端受集中力p=200，左端固定，不考慮梁自重，彈性模量E=210000 Pa，剪切模量v=0.25。采用解析法、自然單元法(采用結(jié)點(diǎn)間距為100的均勻和非均勻兩種布點(diǎn)方案離散懸臂梁，如圖7、8)、有限元法(采用單元尺寸為100的四面體單元和六面體單元來(lái)離散懸臂梁，如圖9、10)等方法計(jì)算梁中心軸線上結(jié)點(diǎn)處的撓度曲線，如圖11所示。

圖6 三維懸臂梁幾何圖形

圖7 三維懸臂梁均勻布點(diǎn)

圖8 三維懸臂梁非均勻布點(diǎn)

圖9 三維懸臂梁六面體單元?jiǎng)澐?/p>

圖10 三維懸臂梁四面體單元?jiǎng)澐?/p>

圖11 懸臂梁沿y=0，z=0軸線上y方向位移比較

從圖11可知，本算例在相同單元尺寸、結(jié)點(diǎn)間距的網(wǎng)格劃分或布點(diǎn)下，四面體單元有限元法的結(jié)點(diǎn)最多，但是撓度值是解析解的70%左右，而均勻布點(diǎn)的自然單元法的計(jì)算精度比四面體精度明顯提高，而且結(jié)點(diǎn)數(shù)也少于四面體單元。自然單元法的均勻布點(diǎn)和非均勻布點(diǎn)的計(jì)算結(jié)果都與解析解非常接近，精度與六面體單元在同一數(shù)量級(jí)上，如表2所示。對(duì)比自然單元法的兩種布點(diǎn)方案，均勻布點(diǎn)比非均勻布點(diǎn)精度高，但是相差不大也在同一數(shù)量級(jí)上。再進(jìn)一步研究自然單元法精度提高的原因發(fā)現(xiàn)，自然單元法的形函數(shù)中對(duì)于任意積分點(diǎn)搜索到自然鄰接點(diǎn)在6個(gè)到20個(gè)之間，最普遍的是8個(gè)鄰接點(diǎn)，并且這些鄰接點(diǎn)多數(shù)都均勻分布在積分點(diǎn)周?chē)?，再根?jù)連續(xù)物質(zhì)的性質(zhì)知，由這些鄰接點(diǎn)插值出的形函數(shù)就比較接近連續(xù)物質(zhì)真實(shí)形函數(shù)。這也就是為什么自然單元法的精度會(huì)與六面體在同一數(shù)量級(jí)原因。如果把四面體單元上進(jìn)行結(jié)點(diǎn)加密，它也可以得到與自然單元法、六面體單元相同數(shù)量級(jí)的精度，見(jiàn)圖12所示。

表2 自然單元法與六面體法數(shù)據(jù)對(duì)較表

圖12 不同插值點(diǎn)下的計(jì)算位移比較

3.3 懸臂梁端部受集中動(dòng)載荷作用

懸臂梁尺寸如下：l=10 m，b=h=1 m，如圖13。左端固定，不考慮梁自重，彈性模量E=210000 Pa，剪切模量v=0.25，密度ρ=2000 kg/m3，阻尼系數(shù)α=0.2、β=0.2。右端突然施加y向集中力p=g(t)，如圖14，沖擊作用時(shí)間t=4 s，時(shí)間步長(zhǎng)Δt=0.01 s。為說(shuō)明本文方法可行性，本文采用彈性小變形下，有限元軟件Ansys的計(jì)算結(jié)果作為精確解，與本文方法相比較說(shuō)明本文正確性。自然單元法(采用結(jié)點(diǎn)間距為0.5的均勻布點(diǎn)，如圖15)、有限元法(采用單元尺寸為0.1六面體單元來(lái)離散懸臂梁，如圖16)等方法分別計(jì)算無(wú)阻尼和有阻尼情況下，梁中心軸線上結(jié)點(diǎn)處的波形曲線，以及軸線中點(diǎn)和自由端點(diǎn)的時(shí)程曲線。

圖13 三維懸臂梁幾何圖形

圖14 輸入荷載圖波形

圖15 三維懸臂梁均勻布點(diǎn)

圖16 三維懸臂梁六面體單元?jiǎng)澐?/p>

由圖17可知，懸臂梁軸線在不同時(shí)間下，計(jì)算得到的無(wú)阻尼波形曲線與有限元法計(jì)算得到曲線非常接近。再由圖18可知，兩種方法計(jì)算的梁端點(diǎn)和中點(diǎn)的位移時(shí)程曲線也非常接近，因此說(shuō)明本文方法在無(wú)阻尼情況下，能夠像有限元法一樣模擬彈性體的動(dòng)力學(xué)問(wèn)題。由圖19可知，在有阻尼條件下，本文方法也可以模擬彈性體動(dòng)力學(xué)問(wèn)題。

圖17 無(wú)阻尼條件下軸線上y方向波形曲線

圖18 無(wú)阻尼條件下點(diǎn)(10,0.5,0.5)y位移時(shí)程曲線

圖19 有阻尼條件下點(diǎn)(10,0.5,0.5)y位移時(shí)程曲線

4 結(jié) 論

本文推導(dǎo)三維自然單元法的動(dòng)力學(xué)問(wèn)題的離散格式，空間域上采用本文推導(dǎo)的三維自然單元法的插值方法進(jìn)行離散，時(shí)間域采用直接積分方法中的中心差分和Newmark常平均加速度法相結(jié)合的第一種積分格式進(jìn)行離散。該積分方法可以消除加速度項(xiàng)，直接利用速度和位移聯(lián)立方程組求解結(jié)點(diǎn)位移和速度。通過(guò)這種解耦方法提高了本文方法的計(jì)算效率。最后以分片試驗(yàn)、懸臂梁為算例，采用大型數(shù)值軟件ANSYS模擬懸臂梁動(dòng)力響應(yīng)的計(jì)算結(jié)果作為基準(zhǔn)，與本文推導(dǎo)的算法相對(duì)比，來(lái)驗(yàn)證本文算法的有效性和合理性。

對(duì)比本文方法和有限元方法，在彈性小變形條件下，兩種方法都可以模擬計(jì)算彈性體動(dòng)力學(xué)問(wèn)，但是在大變形條件下，有限元方法將出現(xiàn)網(wǎng)格畸變使計(jì)算無(wú)法進(jìn)行，然而本文方法則不會(huì)出現(xiàn)網(wǎng)格畸變，計(jì)算仍然可以進(jìn)行。通過(guò)對(duì)比兩種方法的計(jì)算結(jié)果可以知道，本文方法在很少結(jié)點(diǎn)情況下，可以得到與有限元數(shù)倍結(jié)點(diǎn)的計(jì)算結(jié)果相近。由此可看出，本方法受結(jié)點(diǎn)離散距離的影響較小，不像有限元在計(jì)算動(dòng)力學(xué)問(wèn)題時(shí)，受單元尺寸的限制，即單元尺寸要足夠小，否則不能很好反映彈性體的動(dòng)力特性。

[1] Dai K Y, Liu G R, Lim K M, et al. A mesh-free method for static and free vibration analysis of shear deformable laminated composite plates[J]. Journal of Sound and Vibration,2004,269(3-5):633-652.

[2] Braun J，Sambridge M. A numerical method for solving partial differential equations on highly irregular evolving grids[J].Nature，1995，376：655-660.

[3] Sukumar N. The natural element method in solid mechanics[D].Theoretical and Applied Mechanics，Northwestern University，1998.

[4] Sukumar N，Moran B，Belytschko T. The natural element method insolid mechanics[J]. International Journal for Numerical Methods in Engineering，1998，43(5)：839-887.

[5] Sukumar N，Moran B A. Semenov Y，et al. Natural neighbour Galerkin methods[J]. International Journal for Numerical Methods in Engineering，2001，50：1-27.

[6] Cueto E，Calvo B，Doblar′e M. Modeling three-dimensional piece-wise homogeneous domains using the α-shape based natural element method[J]. International Journal for Numerical Methods in Engineering，2002，54(6)：871-897.

[7] 蔡永昌，朱合華，王建華. 基于Voronoi結(jié)構(gòu)的無(wú)網(wǎng)格局部Petrov-Galerkin方法[J].力學(xué)學(xué)報(bào)，2003，35(2)：187-193.

[8] 盧 波，葛修潤(rùn)，王水林.自然單元法數(shù)值積分方案研究[J] 巖石力學(xué)與工程學(xué)報(bào)， 2005，24(11)：1917-1924.

[9] 盧 波，葛修潤(rùn)，王水林.自適應(yīng)自然單元法研究—誤差估計(jì)[J].巖石力學(xué)與工程學(xué)報(bào)，2005，24(22)：4065-4072.

[10]Martniéz M A，Cueto E，Doblaré M，et al. A meshless simulation of injection processes involving short fibers molten composites[J].International Journal of Forming Processes，2001，4(3)：217-236.

[11]朱合華，楊寶紅，蔡永昌，等.無(wú)網(wǎng)格自然單元法在彈塑性分析中的應(yīng)用[J].巖土力學(xué)，2004，25(4)：671-674.