房霆宸， 李 武

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

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

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

1 位移函數(shù)

1.1 近似位移函數(shù)

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

圖1 8結點一階Voronoi圖

標號為1～8的8個結點構成了待插值點X的自然相鄰結點，待插值點與自然相鄰點一起形成Voronoi圖稱為待插值點的二次Voronoi結構。因為空間結構復雜，所以采用1/8結構進行說明推導過程。

圖2 8結點二階Voronoi圖

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

(1)

(2)

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

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

(3)

(4)

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

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

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

0≤ΦI(x)≤1

(5)

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

ΦI(xJ)=δIJ

(6)

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

(7)

(8)

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

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

根據(jù)插值函數(shù)式(3)、式(4)求導相應導數(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)是任意變化的空間位移增量。

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

(20)

(21)

(22)

式中，Φ(x)為結點插值的形函數(shù)，x為結點空間幾何坐標。

將空間離散后的位移表達式(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)分別為質量矩陣、阻尼矩陣、剛度矩陣及荷載向量。

3 算例分析

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

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

Lu=

(28)

3.1 分片試驗

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

圖3 立方體

圖4 均勻布點

圖5 非均勻布點

編號a位移值x方向a坐標值x方向b位移值x方向b坐標值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可知本文三維自然單元法的算法可以精確通過分片試驗，因此該算法解的收斂性通過分片試驗的檢驗，說明本文算法是正確的、可行的。

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。采用解析法、自然單元法(采用結點間距為100的均勻和非均勻兩種布點方案離散懸臂梁，如圖7、8)、有限元法(采用單元尺寸為100的四面體單元和六面體單元來離散懸臂梁，如圖9、10)等方法計算梁中心軸線上結點處的撓度曲線，如圖11所示。

圖6 三維懸臂梁幾何圖形

圖7 三維懸臂梁均勻布點

圖8 三維懸臂梁非均勻布點

圖9 三維懸臂梁六面體單元劃分

圖10 三維懸臂梁四面體單元劃分

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

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

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

圖12 不同插值點下的計算位移比較

3.3 懸臂梁端部受集中動載荷作用

懸臂梁尺寸如下：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，沖擊作用時間t=4 s，時間步長Δt=0.01 s。為說明本文方法可行性，本文采用彈性小變形下，有限元軟件Ansys的計算結果作為精確解，與本文方法相比較說明本文正確性。自然單元法(采用結點間距為0.5的均勻布點，如圖15)、有限元法(采用單元尺寸為0.1六面體單元來離散懸臂梁，如圖16)等方法分別計算無阻尼和有阻尼情況下，梁中心軸線上結點處的波形曲線，以及軸線中點和自由端點的時程曲線。

圖13 三維懸臂梁幾何圖形

圖14 輸入荷載圖波形

圖15 三維懸臂梁均勻布點

圖16 三維懸臂梁六面體單元劃分

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

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

圖18 無阻尼條件下點(10,0.5,0.5)y位移時程曲線

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

4 結 論

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

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

[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結構的無網(wǎng)格局部Petrov-Galerkin方法[J].力學學報，2003，35(2)：187-193.

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

[9] 盧 波，葛修潤，王水林.自適應自然單元法研究—誤差估計[J].巖石力學與工程學報，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ǎng)格自然單元法在彈塑性分析中的應用[J].巖土力學，2004，25(4)：671-674.