基于點基局部光滑點插值法的結構動力分析

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

(1.大連理工大學 船舶工程學院 遼寧省深海浮動結構工程實驗室，大連 116024；2.大連理工大學工業裝備與結構分析國家重點實驗室，大連 116024；3.高新船舶與深海開發裝備協同創新中心，上海 200240；4.大連海洋大學 航海與船舶工程學院，大連 116023)

經過半個世紀的發展，有限元法已成為科學研究以及工程應用方面強有力的數值工具，并擁有諸多的商業軟件，能夠對各類復雜物理問題進行分析求解，但其同時也存在一些固有缺陷：其中之一就是三角形單元(四面體)網格劃分簡便，但其剛度過硬，導致計算精度較差和收斂性不佳；四邊形(六面體)單元精度高但難以實現復雜構件的網格自動剖分。

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

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

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

1 數值模型

1.1 有限元法

傳統有限元基于如下伽遼金弱形式

(1)

位移場被近似為

(2)

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

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

KFEMd=f

(3)

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

(4)

(5)

式中

(6)

本文采用三節點三角形單元，應變矩陣是常數矩陣，則式(4)變為

(7)

式中:Ve為單元面積。

1.2 點基光滑點插值法

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

圖1 三角形背景網格以及在此基礎上形成的基于節點的應變光滑域

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

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

(8)

(9)

(10)

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

位移場近似表示為

(11)

(12)

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

(13)

(14)

(15)

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

(16)

式中

(17)

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

(18)

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

(19)

式中

(20)

(21)

(22)

(23)

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

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

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

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

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

(24)

式中:ρ為質量密度。

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

(25)

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

(26)

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

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

(27)

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

(28)

(29)

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

C=αM+βK

(30)

式中:α,β為阻尼系數。

(31)

(32)

(33)

1.4 數值實現

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

1) 對節點光滑域循環。

2) 對光滑域中的單元循環。

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

(a) 根據節點選擇方案選擇插值節點并計算點插值函數。

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

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

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

4) 結束高斯點循環。

5) 結束單元循環。

6) 結束節點光滑域循環。

7) 對單元循環。

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

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

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

8) 結束單元循環。

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

10) 施加本質邊界條件。

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

2 分片試驗

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

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

(34)

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

(35)

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

3 數值算例

3.1 細長懸臂板的自由振動

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

表1列舉了采用相同節點離散、不同方法求得的該懸臂板的前八階固有頻率值，圖4給出了前八階模態。根據數據結果可知，當前的NPS-PIM方法求得的

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

Tab.1Firsteightnaturalfrequencies(in10kHz)oftheslendercantilever

NS-PIM NPS-PIM 有限元(T3) 參考解(FEM-Q4)網格:448 0.080 41 0.083 85 0.085 35 0.082 4節點: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網格:1 074 0.081 47 0.082 86 0.083 50 0.082 4節點: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網格:2 004 0.081 83 0.082 56 0.082 91 0.082 4 節點: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

模態1，f=838.5 Hz

模態2，f=5 022 Hz

模態3,f=12 830 Hz

模態4，f=13 210 Hz

模態5，f=23 950 Hz

模態6，f=36 460 Hz

模態7，f=38 430 Hz

模態8，f=50 100 Hz

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

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

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

3.2 剪切墻的自由振動

如圖5所示，該算例將研究有四個開口的剪切墻，墻底部完全固定，計算狀態為平面應力。利用522個不規則節點離散該問題域。基本數據為E=1.0×104Pa，v=0.2，t=1.0 m，ρ=1.0 kg/m3。

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

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

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

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

圖6 NPS-PIM求得的剪切墻的前八階模態

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

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

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

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

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

圖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點出的瞬態位移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 結 論

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

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

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

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

(4)在采用同樣線性三角形單元網格對問題域進行離散的情況下，NPS-PIM在求解自由振動以及受迫振動時均較有限元準確，基于能夠實現自動剖分的三角形背景網格亦能達到較高精度。