張桂勇， 魯 歡， 王海英， 于大鵬， 孫鐵志

(1.大連理工大學(xué) 船舶工程學(xué)院 遼寧省深海浮動結(jié)構(gòu)工程實驗室，大連 116024；2.大連理工大學(xué)工業(yè)裝備與結(jié)構(gòu)分析國家重點實驗室，大連 116024；3.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240；4.大連海洋大學(xué) 航海與船舶工程學(xué)院，大連 116023)

經(jīng)過半個世紀(jì)的發(fā)展，有限元法已成為科學(xué)研究以及工程應(yīng)用方面強(qiáng)有力的數(shù)值工具，并擁有諸多的商業(yè)軟件，能夠?qū)Ω黝悘?fù)雜物理問題進(jìn)行分析求解，但其同時也存在一些固有缺陷：其中之一就是三角形單元(四面體)網(wǎng)格劃分簡便，但其剛度過硬，導(dǎo)致計算精度較差和收斂性不佳；四邊形(六面體)單元精度高但難以實現(xiàn)復(fù)雜構(gòu)件的網(wǎng)格自動剖分。

為了克服有限元法的缺陷，國內(nèi)外學(xué)者進(jìn)行了大量的研究。Pian等[1]提出了混合有限元法以改善傳統(tǒng)有限元法的計算性能；Liu等[2-4]在提出G空間理論和廣義梯度光滑技術(shù)的基礎(chǔ)上，進(jìn)一步提出了無網(wǎng)格光滑點插值法(S-PIM)，點基光滑點插值法(NS-PIM)就是其中的一種。研究發(fā)現(xiàn)點基光滑點插值法具有諸多良好的性質(zhì)：應(yīng)力/應(yīng)變解準(zhǔn)確，避免剪切鎖死，更為重要的是其能給出能量解的上界。已知有限元結(jié)果能夠提供能量解下界，通過點基光滑點插值法獲得能量解上界，從而允許我們從數(shù)值角度獲得復(fù)雜問題能量范數(shù)的精確解區(qū)間[5-6]。

但點基光滑點插值法有其自身缺陷，表現(xiàn)為剛度“過軟”[7-8]。研究發(fā)現(xiàn)：其在動力學(xué)求解時，常會出現(xiàn)虛假的非零能模式；盡管采用無條件收斂的直接積分法如Newmark法[9]，時間也會不穩(wěn)定[10-11]；且在求解固有頻率時，常常會出現(xiàn)計算結(jié)果偏低，在高階頻率時表現(xiàn)尤為明顯。

針對點基光滑點插值法的上述缺陷，本文將其與有限元法結(jié)合，對問題域進(jìn)行局部應(yīng)變光滑，借助著有限元法剛度“過硬”特性[12]，來矯正點基光滑點插值法剛度“過軟”情況，以期得到接近實際剛度的模型，并改善原始點基光滑點插值法動力求解時的時間不穩(wěn)定狀況。為了驗證該方法的有效性，進(jìn)行了一系列固有頻率和響應(yīng)計算。結(jié)果表明本文所提出方法能夠有效克服點基光滑點插值法剛度過軟造成的虛假非零能模式和固有頻率過低的問題，在采用同樣線性三角形單元網(wǎng)格對問題域進(jìn)行離散的情況下，能夠給出較有限元模型更為準(zhǔn)確的數(shù)值結(jié)果。

1 數(shù)值模型

1.1 有限元法

傳統(tǒng)有限元基于如下伽遼金弱形式

(1)

位移場被近似為

(2)

式中:NJ(x)為形函數(shù)矩陣;dJ為節(jié)點位移矢量;NP為插值所用的節(jié)點數(shù)。

將式(2)代入式(1)中得到系統(tǒng)的離散方程

KFEMd=f

(3)

式中:KFEM為有限元系統(tǒng)剛陣;f為力矢量，由以下公式得到

(4)

(5)

式中

(6)

本文采用三節(jié)點三角形單元，應(yīng)變矩陣是常數(shù)矩陣，則式(4)變?yōu)?/p>

(7)

式中:Ve為單元面積。

1.2 點基光滑點插值法

點基光滑點插值法借助背景單元對問題域離散，在此基礎(chǔ)上形成基于節(jié)點的應(yīng)變光滑域，如圖1所示，Ω(K)為節(jié)點K所在的光滑域，Γ(K)為光滑域邊界。

圖1 三角形背景網(wǎng)格以及在此基礎(chǔ)上形成的基于節(jié)點的應(yīng)變光滑域

Fig.1 Triangular background elements and strain smoothing domains associated with field nodes

根據(jù)GS-Galerkin(Generalized smoothed Galerkin weak form)弱式得到下式

(8)

(9)

(10)

式中:n1和n2分別對應(yīng)光滑域邊界的方向余弦。

位移場近似表示為

(11)

(12)

將式(10)、(11)代入式(9)中得到

(13)

(14)

(15)

將式(10)和形函數(shù)ΦI(x)代入式(15)中得到

(16)

式中

(17)

式中:nl(x)為Ln矩陣中的元素。采用高斯積分對式(17)計算，得到

(18)

將式(11)、(13)代入式(8)中得到

(19)

式中

(20)

(21)

(22)

(23)

1.3 點基局部光滑點插值法求解動力問題

由前面介紹得知，點基光滑點插值法將每個三角形背景網(wǎng)格分為三部分，參與以對應(yīng)節(jié)點為中心的梯度光滑過程，并且前面的研究已經(jīng)證明此方法剛度過軟；而有限元法則可看作每個背景單元均未進(jìn)行梯度光滑過程，其剛度過硬。綜合以上兩種方法的特點，如圖2所示，本文方法把背景單元的每條邊三等分，每個等分點分別連接相鄰兩邊最近的等分點，將背景單元分割為頂點處三個三角形區(qū)域和中間空白六邊形區(qū)域。其中三個三角形區(qū)域(圖2中陰影部分)參與形成對應(yīng)節(jié)點光滑域，而中間的六邊形區(qū)域則不進(jìn)行梯度光滑，因此方法可以看作是基于點的局部梯度光滑點插值法。

圖2 點基局部光滑點插值法問題域背景網(wǎng)格及光滑域形成

Fig.2 Background elements and smoothing domains of the NPS-PIM

對于動力問題GS-Galerkin弱式可寫為

(24)

式中:ρ為質(zhì)量密度。

將各個物理量代入式(24)得到

(25)

此式為動力學(xué)求解的一般形式，其中，M為質(zhì)量矩陣，C為阻尼矩陣，f為力矢量，因此，只要得到各個矩陣，求解便可以得到固有頻率值和響應(yīng)值，本文點基局部光滑點插值法的剛度陣是將有限元法和點基光滑點插值法的二者的剛度陣結(jié)合，剛度陣K的求解如下

(26)

質(zhì)量陣M采用的是集中質(zhì)量陣，是將每個單元的質(zhì)量平均地分配給該單元的節(jié)點，對所有單元循環(huán)組裝得到總的質(zhì)量陣。

對于自由振動，不考慮外力和阻尼，因此式(25)簡化為

(27)

求解該式有時需要施加邊界條件，大多為本質(zhì)邊界，本文計算采用消去法，直接將邊界條件節(jié)點對應(yīng)的剛度陣和質(zhì)量陣處的元素劃去。式(27)的解可表示為

(28)

(29)

對于受迫振動，需要考慮阻尼的影響并施加外力，本文中的阻尼采用瑞利阻尼

C=αM+βK

(30)

式中:α,β為阻尼系數(shù)。

(31)

(32)

(33)

1.4 數(shù)值實現(xiàn)

點基局部光滑點插值法求解自由振動問題時數(shù)值流程如下：

1) 對節(jié)點光滑域循環(huán)。

2) 對光滑域中的單元循環(huán)。

3) 對單元中光滑域邊的高斯點循環(huán)。

(a) 根據(jù)節(jié)點選擇方案選擇插值節(jié)點并計算點插值函數(shù)。

(b) 計算應(yīng)變位移矩陣。

(c) 計算節(jié)點剛度陣及力矢量。

(d) 將節(jié)點剛度陣及力矢量組裝到整體剛陣及整體力矢量。

4) 結(jié)束高斯點循環(huán)。

5) 結(jié)束單元循環(huán)。

6) 結(jié)束節(jié)點光滑域循環(huán)。

7) 對單元循環(huán)。

(a) 計算單元梯度矩陣B。

(b) 計算單元剛度陣、質(zhì)量陣以及力矢量。

(c) 將單元剛度陣、質(zhì)量陣以及力矢量組裝。

8) 結(jié)束單元循環(huán)。

9) 將有限元和點基光滑點插值法的總剛疊加，得到最終的剛度矩陣。

10) 施加本質(zhì)邊界條件。

11) 根據(jù)最終剛度陣和質(zhì)量陣求解得到各階固有頻率。

2 分片試驗

檢驗一個數(shù)值方法是否收斂于真實解可利用分片試驗，其要求內(nèi)部與邊界的節(jié)點位移滿足相同的線性函數(shù)[13](機(jī)器誤差范圍內(nèi))。

如圖3所示，這里利用132個節(jié)點離散該方片，以探究本文方法是否收斂于真實解，加載邊界處的線性位移為

(34)

上述位移亦為方片的解析解。利用位移誤差范數(shù)來檢驗位移誤差，如下式

(35)

式中，上標(biāo)“nume” 表示數(shù)值解，上標(biāo)“exact”表示解析解，Nn為總的節(jié)點數(shù)。計算得到本文方法的上述位移誤差范數(shù)為0.122 5×10-14，證明了該方法能夠通過分片試驗。

3 數(shù)值算例

3.1 細(xì)長懸臂板的自由振動

本節(jié)將研究細(xì)長懸臂板的自由振動，其尺寸為：長L=100 m，寬D=10 m，厚度t=1.0 m，楊氏模量2.1×104Pa，泊松比v=0.3，質(zhì)量密度ρ=8.0×10-10kg/m3。利用三種不同的網(wǎng)格劃分離散該問題域，分別為：274，612，和1 104個不規(guī)則節(jié)點。由于歐拉理論懸臂梁解析解未考慮剪切變形的影響，所以本文利用精細(xì)網(wǎng)格劃分100×10的FEM-Q4單元求解該問題得到的解作為參考解[14]。

表1列舉了采用相同節(jié)點離散、不同方法求得的該懸臂板的前八階固有頻率值，圖4給出了前八階模態(tài)。根據(jù)數(shù)據(jù)結(jié)果可知，當(dāng)前的NPS-PIM方法求得的

表1細(xì)長懸臂梁前八階固有頻率(10kHz)

Tab.1Firsteightnaturalfrequencies(in10kHz)oftheslendercantilever

NS-PIM NPS-PIM 有限元(T3) 參考解(FEM-Q4)網(wǎng)格:448 0.080 41 0.083 85 0.085 35 0.082 4節(jié)點:274 0.481 2 0.502 2 0.511 0 0.494 4 1.264 1.283 1.283 1.282 4 1.282 1.321 1.343 1.302 2 2.288 2.395 2.434 2.366 3 3.476 3.646 3.705 3.608 5 3.841 3.843 3.844 3.844 2 4.748 5.010 5.088 4.967 4網(wǎng)格:1 074 0.081 47 0.082 86 0.083 50 0.082 4節(jié)點:612 0.488 3 0.496 7 0.500 5 0.494 4 1.282 1.282 1.283 1.282 4 1.285 1.307 1.317 1.302 2 2.332 2.373 2.391 2.366 3 2.492 3.616 3.643 3.608 5 3.246 3.844 3.844 3.844 2 3.550 4.975 5.012 4.967 4網(wǎng)格:2 004 0.081 83 0.082 56 0.082 91 0.082 4 節(jié)點:1 104 0.490 6 0.495 1 0.497 2 0.494 4 1.2821.2821.2821.282 4 1.2911.3031.3091.302 2 2.3452.3672.3772.366 3 3.5743.6093.6233.608 5 3.8443.8443.8443.844 2 4.9164.9664.9874.967 4

模態(tài)1，f=838.5 Hz

模態(tài)2，f=5 022 Hz

模態(tài)3,f=12 830 Hz

模態(tài)4，f=13 210 Hz

模態(tài)5，f=23 950 Hz

模態(tài)6，f=36 460 Hz

模態(tài)7，f=38 430 Hz

模態(tài)8，f=50 100 Hz

圖4 NPS-PIM方法求得的懸臂板前八階振動模態(tài)

Fig.4 First 8 free vibration modes of the slender cantilever by NPS-PIM

固有頻率比有限元小，并克服了點基光滑點插值法解固有頻率值過低的缺陷，較之二者均更接近真實解，且前八階模態(tài)中沒有虛假的非零能模態(tài)出現(xiàn)，時間穩(wěn)定。

3.2 剪切墻的自由振動

如圖5所示，該算例將研究有四個開口的剪切墻，墻底部完全固定，計算狀態(tài)為平面應(yīng)力。利用522個不規(guī)則節(jié)點離散該問題域?；緮?shù)據(jù)為E=1.0×104Pa，v=0.2，t=1.0 m，ρ=1.0 kg/m3。

圖5 帶有四個開口的剪切墻及背景網(wǎng)格劃分Fig.5 A shear wall with four openings and its background meshes

在相同節(jié)點離散下、用不同方法求得的前八階固有頻率列于表2中，相應(yīng)地用本文方法求得的模態(tài)繪于圖6中。從結(jié)果中可知，在相同的節(jié)點離散下，較之有限元法(T3)，本文提出的NPS-PIM方法其結(jié)果和參考解更加接近；與點基光滑點插值方法相比，克服了其求解固有頻率值過低的缺陷，更接近參考解。并由圖6可以看出NPS-PIM計算得到的前八階陣型并未出現(xiàn)虛假的非零能模態(tài)，克服了原始點基光滑點插值法的時間不穩(wěn)定性。

表2 剪切墻的前八階固有頻率Tab.2 First eight frequencies(rad/s) of a shear wall rad/s

模態(tài)1f=2.075 rad/s模態(tài)2f=7.101 rad/s模態(tài)3f=7.625 rad/s模態(tài)4f=11.95 rad/s模態(tài)5f=15.36 rad/s模態(tài)6f=18.34 rad/s模態(tài)7f=19.88 rad/s模態(tài)8f=22.21 rad/s

圖6 NPS-PIM求得的剪切墻的前八階模態(tài)

Fig.6 First 8 free vibration modes of the shear wall by the NPS-PIM

3.3 端部施加簡諧載荷的細(xì)長懸臂板的受迫振動

如圖7所示，受迫振動的算例是一個懸臂板，計算狀態(tài)為平面應(yīng)力狀態(tài)。其基本數(shù)據(jù)為：楊氏模量E=3×107Pa，泊松比v=0.3，厚度t=1.0 m，時間間隔Δt=1×10-3s。該算例中，懸臂板端部受到一個簡諧載荷P=1 000g(t)。計算時，利用275個不規(guī)則節(jié)點離散該問題域。為簡單起見，取質(zhì)量密度ρ=1.0 kg/m3，取懸臂梁端部點A的Y方向位移作為考察。為了比較，利用有限元FEM-T3在相同節(jié)點離散下計算該問題，并以有限元FEM-Q4精細(xì)網(wǎng)格結(jié)果作為參考解。

首先，考慮載荷g(t)=sin(ωft)的情況，其中，ωf=27是載荷的頻率。本文采用瑞利阻尼，其兩個參數(shù)分別為α，β。對于振動方程的求解，這里采用Newmark法，當(dāng)0.5≤θ≤1，Newmark法無條件穩(wěn)定，現(xiàn)取θ=0.5。

如圖8所示，采用相對較大的時間步長Δt=5×10-3s，計算20 s內(nèi)的位移響應(yīng)，可以看出，NPS-PIM方法求得的位移響應(yīng)十分穩(wěn)定。如圖9所示，采用較小的時間步長Δt=1×10-3s，分別計算時程0.2 s和2 s時的響應(yīng)，從圖中可以看出，本文NPS-PIM法比有限元(T3)計算的響應(yīng)結(jié)果更為準(zhǔn)確，更接近參考解(有限元Q4)。如圖10所示，在無阻尼(α=0,β=0)，懸臂板受到大小為g(t)=1持續(xù)時間0.5 s的瞬態(tài)載荷時，分別計算了時程為2 s和20 s的響應(yīng)，從圖中可以看出，NPS-PIM求得的瞬態(tài)振動的解較有限元精確，在有阻尼的情況下(α=0.000 5，β=0.000 7)求得的響應(yīng)亦十分穩(wěn)定。

圖8 NPS-PIM求得的A點處的位移uy(g(t)=sin(ωt))Fig.8 Displacement uy at the point A using NPS-PIM(g(t)=sin(ωt))(a) 時間歷程0.2 s(b) 時間歷程2 s圖9 不同方法求得的點A處的位移uy(g(t)=sin(ωt))Fig.9 Displacements uy at the point A using different methods (g(t)=sin(ωt))

(a) 無阻尼

(b) 在有阻尼的情況下圖10 不同的方法求得的A點出的瞬態(tài)位移uy(α=0.000 5，β=0.000 7)Fig.10 Transient displacements uy at the point A using different methods (α=0.000 5，β=0.000 7)

4 結(jié) 論

本文將有限元法和點基光滑點插值法結(jié)合(NPS-PIM)，充分利用有限元法過剛、點基光滑點插值法過軟的特性，得到了更接近實際剛度的計算模型，并利用數(shù)值算例驗證了該方法的準(zhǔn)確性。通過這些算例我們可以得到以下結(jié)論：

(1)NPS-PIM方法能夠通過分片試驗，收斂于真實解。

(2)NPS-PIM方法克服了NS-PIM的時間不穩(wěn)定性，求解動力問題時沒有出現(xiàn)虛假的非零能模態(tài)，時間穩(wěn)定。

(3)NPS-PIM方法克服了點基光滑點插值法求解固有頻率值過低的缺陷，得到的固有頻率值更為接近參考解。

(4)在采用同樣線性三角形單元網(wǎng)格對問題域進(jìn)行離散的情況下，NPS-PIM在求解自由振動以及受迫振動時均較有限元準(zhǔn)確，基于能夠?qū)崿F(xiàn)自動剖分的三角形背景網(wǎng)格亦能達(dá)到較高精度。