結(jié)合虛擬點(diǎn)的單元基光滑點(diǎn)插值法在結(jié)構(gòu)動力學(xué)分析中的應(yīng)用

張桂勇,陶東松,陳澤聰,陳毓珍

(1.大連理工大學(xué)船舶工程學(xué)院,116024,遼寧大連;2.大連理工大學(xué)工業(yè)裝備與結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室,116024,遼寧大連;3.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心,200240,上海;4.中國船舶工業(yè)集團(tuán)上海船舶研究設(shè)計(jì)院,201203,上海)

在工程問題中,結(jié)構(gòu)的振動固有頻率和動態(tài)響應(yīng)是結(jié)構(gòu)分析的基礎(chǔ),受到了人們的廣泛重視。隨著有限元法[1]在結(jié)構(gòu)數(shù)值模擬方面的迅速發(fā)展,多種商業(yè)數(shù)值仿真軟件被開發(fā)出來,成為工程技術(shù)人員進(jìn)行結(jié)構(gòu)分析的有力工具。但是,有限元法也存在一些固有問題[1],比如采用劃分簡便的線性三角形(四面體)單元求解時(shí),得到的系統(tǒng)剛度過硬,精度較差,而對于精度較高的四邊形(六面體)計(jì)算網(wǎng)格,在復(fù)雜的幾何條件下常常難以自動劃分質(zhì)量較高的網(wǎng)格。

為解決上述問題,近年來發(fā)展出了一系列新的數(shù)值方法,如光滑有限元方法(S-FEM)[2-3]和光滑點(diǎn)插值方法(S-PIM)[4]。S-PIM基于G空間理論[5]和廣義梯度光滑技術(shù)(GSM)[6],相較于剛度過硬的有限元方法,可以有效軟化系統(tǒng)的模型剛度[7-10],進(jìn)而提高計(jì)算性能。在各種S-PIM中,單元基光滑點(diǎn)插值法(CS-PIM)[11]的剛度軟化程度介于點(diǎn)基光滑點(diǎn)插值法(NS-PIM)[12]與邊基光滑點(diǎn)插值法(ES-PIM)之間[13],可在僅使用線性三角形(四面體)單元的情況下,得到令人滿意的結(jié)果。

在各種CS-PIM中,采用四點(diǎn)插值模式(T4)的徑向基光滑點(diǎn)插值法(CSRPIM-T4)計(jì)算精度較低。為提高計(jì)算精度,同時(shí)不增加計(jì)算成本,將一種引入虛擬點(diǎn)參與構(gòu)造的濃縮RPIM形函數(shù)(Cd)[14]應(yīng)用于上述方法,得到結(jié)合虛擬點(diǎn)的單元基光滑點(diǎn)插值法(CSRPIM-T4-Cd)。研究表明:在傳統(tǒng)的虛擬點(diǎn)布置方式下系統(tǒng)剛度過軟[14],導(dǎo)致求得的固有頻率偏低并出現(xiàn)虛假的非零能模態(tài);在動力學(xué)求解中,即使采用無條件收斂的直接積分法,如Newmark法,依然存在時(shí)間不穩(wěn)定性問題。

針對傳統(tǒng)虛擬點(diǎn)布置中存在的問題,本文通過改變虛擬點(diǎn)的布置,可使原本過軟的系統(tǒng)剛度更加接近真實(shí)的剛度,從而提高計(jì)算精度。通過對一系列自由振動與受迫振動的算例進(jìn)行計(jì)算,驗(yàn)證了該方法的精確性。

1 單元基光滑點(diǎn)插值法

圖1 基于三角形網(wǎng)格的單元基光滑域

根據(jù)GS-Galerkin弱形式(generalized smoothed Galerkin weak form)[4],可得

(1)

(2)

(3)

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

位移場表示為

(4)

式中:Sn為用于插值的節(jié)點(diǎn)集;φI(x)為節(jié)點(diǎn)I處的形函數(shù);Φ為由φI(x)組成的形函數(shù)矩陣;dI為節(jié)點(diǎn)I處的位移;d為由dI組成的節(jié)點(diǎn)位移矢量。

將式(3)、式(4)代入式(2),得

(5)

(6)

其中的各元素通過下式得到

(7)

采用高斯積分對式(7)進(jìn)行計(jì)算

(8)

將式(4)和式(5)代入式(1),得

(9)

式中

(10)

(11)

(12)

式中Ne為背景網(wǎng)格的單元數(shù)目。因?yàn)楸尘熬W(wǎng)格同時(shí)作為局部應(yīng)變光滑域,所以光滑域數(shù)目等于背景網(wǎng)格數(shù)。

2 濃縮形函數(shù)

2.1 RPIM形函數(shù)

采用添加了多項(xiàng)式基的徑向基函數(shù)(RBF)作為插值形式,在局部區(qū)域內(nèi)選取支持點(diǎn)進(jìn)行插值模擬。此插值方法可使徑向基PIM(RPIM)具備線性再生性,提高插值精度,同時(shí)避免了多項(xiàng)式基PIM的奇異性問題,方便支持點(diǎn)的靈活選取。使用如下插值形式來模擬位移場[15]

RT(x)a+PT(x)b

(13)

式中:n為局部區(qū)域內(nèi)的支持點(diǎn)數(shù);m為多項(xiàng)式基的項(xiàng)數(shù),此處取線性完全多項(xiàng)式,即m=3;Ri(x)和ai分別表示徑向基及其系數(shù),Pj(x)和bj分別表示多項(xiàng)式基及其系數(shù),形式如下

(14)

(15)

其中αc、q為復(fù)合二次形式徑向基函數(shù)(MQ-RBF)的形狀參數(shù),取文獻(xiàn)[16]中的推薦值,ri為支持點(diǎn)i和被插值點(diǎn)的徑向距離,dc為問題域內(nèi)節(jié)點(diǎn)的平均間距。

式(13)中有n+m個變量,若要使用支持域中n個點(diǎn)的值求解該式,需添加下面m個約束條件作為附加方程

PT(x)a=0

(16)

最后,求解得到對應(yīng)于局部區(qū)域內(nèi)n個節(jié)點(diǎn)位移的RPIM形函數(shù)[15],形式如下

Φ(x)={φ1(x),φ2(x),…,φn(x)}

(17)

2.2 支持點(diǎn)的選取方式

2.2.1 邊基T4選點(diǎn)模式 在單元基光滑點(diǎn)插值法中,光滑域?yàn)槿切伪尘皢卧?光滑域的邊界即三角形的邊。邊基T4選點(diǎn)模式采用基于邊與場節(jié)點(diǎn)位置關(guān)系的方法來尋找支持點(diǎn)。

如圖2所示,對內(nèi)部邊上的點(diǎn),選取與該邊相鄰的2個三角形的4個頂點(diǎn)作為支持點(diǎn),對問題域邊界邊上的點(diǎn),選取該邊的2個頂點(diǎn)作為支持點(diǎn),可以確保邊界條件的正確施加。

(a)內(nèi)部邊的T4選點(diǎn)模式

(b)邊界邊的T4選點(diǎn)模式圖2 基于邊的T4選點(diǎn)模式

2.2.2 加入虛擬點(diǎn)的邊基T4選點(diǎn)模式 RPIM形函數(shù)常需要較多支持點(diǎn)來保證其精確性[14],但這又會降低計(jì)算效率,而使用虛擬點(diǎn)可以在不增加自由度數(shù)量的前提下提高形函數(shù)的模擬精度。

傳統(tǒng)的虛擬點(diǎn)布置方式為:基于邊基T4選點(diǎn)模式,在相鄰的4條邊的中心引入4個虛擬節(jié)點(diǎn),如圖3a所示,其坐標(biāo)值由下式求得

(18)

(a)加入4個虛擬點(diǎn)

(b)添加6個虛擬點(diǎn)圖3 加入虛擬點(diǎn)的內(nèi)部邊T4選點(diǎn)模式

研究[14]發(fā)現(xiàn),引入4個虛擬點(diǎn)計(jì)算得出的系統(tǒng)剛度陣過軟。本文采用一種新的虛擬點(diǎn)布置方式來改善形函數(shù)的模擬精度,如圖3b所示,基于傳統(tǒng)的虛擬點(diǎn)布置方式,在相鄰的2個單元的形心處再引入2個虛擬點(diǎn)(即點(diǎn)9和10),其坐標(biāo)值為

(19)

2.3 濃縮的RPIM形函數(shù)

2.3.1 CSRPIM-T4-Cd4形函數(shù) 使用4個真實(shí)場節(jié)點(diǎn)(1～4)和4個虛擬點(diǎn)(5～8)共8個支持點(diǎn)進(jìn)行插值,可得濃縮前的RPIM形函數(shù)

Φs(x)={φ1(x),φ2(x),…,φ8(x)}

(20)

計(jì)算點(diǎn)處的位移可以表示為

(21)

虛擬點(diǎn)的位移可以通過真實(shí)場節(jié)點(diǎn)的位移線性插值得到

(22)

將式(22)代入式(21),可以得到

(23)

濃縮后的RPIM形函數(shù)可以表示為如下形式

(24)

2.3.2 CSRPIM-T4-Cd6形函數(shù) 與前面相似,使用4個真實(shí)場節(jié)點(diǎn)(1～4)和6個虛擬點(diǎn)(5～10)共10個支持點(diǎn)進(jìn)行插值,可得濃縮前的RPIM形函數(shù)

Φs(x)={φ1(x),φ2(x),…,φ10(x)}

(25)

虛擬點(diǎn)位移由真實(shí)場節(jié)點(diǎn)位移線性插值得到

(26)

同樣地,可以得到在引入6個虛擬點(diǎn)后的濃縮RPIM形函數(shù)

(27)

3 結(jié)構(gòu)動力問題求解

3.1 系統(tǒng)離散方程

離散的動力問題GS-Galerkin弱形式[4]為

(28)

式中:Ω為問題域;D為材料彈性矩陣;tΓ為自然邊界Γt上的外力向量;b為體力向量;c為阻尼系數(shù);ρ為質(zhì)量密度。

將式(4)和式(5)代入式(28),得

(29)

式中:C為阻尼矩陣;M為質(zhì)量矩陣;f為力矢量。

(30)

(31)

(32)

質(zhì)量矩陣采用動力問題中常用的集中質(zhì)量陣,將單元質(zhì)量集中到各個節(jié)點(diǎn)形成對角陣的形式,使方程便于求解。

3.2 自由振動

對于自由振動,不考慮阻尼和外力的影響,式(29)可簡化為

(33)

式(33)的通解可以寫成

d=dAexp(iωt)

(34)

式中:dA為位移幅值;t為時(shí)間;ω為自然頻率,可由下式求得

(35)

3.3 受迫振動

受迫振動需要考慮阻尼以及外力的影響。采用瑞利阻尼,表達(dá)式為

(36)

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

(37)

(38)

(39)

以上基于二維三角形單元的CSRPIM方法,可以擴(kuò)展到基于四面體單元離散的三維問題[4]。在三維問題中,通過基于四面體單元的梯度光滑技術(shù),將基于四面體的局部積分轉(zhuǎn)化為基于四面體各個三角形面的積分;對于各個積分點(diǎn)的插值,同樣采用徑向基單元點(diǎn)插值法,在基于面的插值點(diǎn)選擇模式基礎(chǔ)上,通過在相鄰四面體邊線及形心布置虛擬點(diǎn),達(dá)到提高插值精度、改變模型剛度、進(jìn)而提高數(shù)值計(jì)算準(zhǔn)確性的目的。

4 分片試驗(yàn)

通過標(biāo)準(zhǔn)分片試驗(yàn)可以驗(yàn)證本文方法能否收斂于準(zhǔn)確解,它要求該分片內(nèi)部的所有節(jié)點(diǎn)在機(jī)器精度上與分片的邊界節(jié)點(diǎn)滿足相同的位移函數(shù)[1]。

使用圖4所示的背景網(wǎng)格來驗(yàn)證本文的數(shù)值方法,單位面積的方片由132個分布節(jié)點(diǎn)劃分成222個單元。

圖4 132個節(jié)點(diǎn)離散的方片

在所有的邊界節(jié)點(diǎn)上施加如下線性位移

ux=0.6x;uy=0.6y

(40)

計(jì)算精度用如下位移誤差模來表示

(41)

式中:上角標(biāo)“ana”表示解析解,“num”表示數(shù)值解;Nn為總節(jié)點(diǎn)數(shù)。

采用標(biāo)準(zhǔn)分片試驗(yàn)計(jì)算得到的位移誤差模如下:CSRPIM-T4-Cd4的位移誤差模為6.834 517 3×10-12,CSRPIM-T4-Cd6的位移誤差模為3.742 340 5×10-11。以上結(jié)果顯示,使用2種濃縮形函數(shù)的CSRPIM均能通過分片試驗(yàn)。

5 數(shù)值算例

5.1 連桿的自由振動

對一連桿進(jìn)行自由振動分析。連桿尺寸如圖5所示,完全約束左側(cè)半徑為1.0 cm的大孔內(nèi)表面。

圖5 連桿的幾何尺寸

采用平面應(yīng)力計(jì)算狀態(tài),材料參數(shù)為:楊氏模量E=10 GPa,泊松比ν=0.3,質(zhì)量密度ρ=7.8×103kg/m3。用含633個不規(guī)則節(jié)點(diǎn)、1 010個單元的三角形背景網(wǎng)格離散該問題域,并應(yīng)用ANSYS有限元軟件以十分細(xì)密的四節(jié)點(diǎn)四邊形單元(21 923個節(jié)點(diǎn),21 146個單元)計(jì)算得到模態(tài)參考解。

表1為3種方法求解出的前10階固有頻率以及各階頻率的相對誤差,同樣可以看出:FEM-T3的固有頻率偏高,系統(tǒng)剛度過硬;CSRPIM-T4-Cd4的固有頻率偏低,系統(tǒng)剛度過軟;CSRPIM-T4-Cd6的精確度最高,最接近真實(shí)剛度。表2是2種虛擬點(diǎn)布置方式求得的固有頻率和模態(tài)振型與參考解的比較,從中可以看到:CSRPIM-T4-Cd4由于剛度過軟,求解得到的固有頻率偏低,且在高階時(shí)存在虛假的非零能模態(tài)(加框振型);CSRPIM-T4-Cd6能夠精確地模擬系統(tǒng)剛度,得到準(zhǔn)確的固有頻率和振型。

表1 3種方法計(jì)算的連桿前10階固有頻率及其誤差

表2 連桿的模態(tài)振型及固有頻率

5.2 懸臂梁的受迫振動

懸臂梁如圖6所示,長L=48 m,寬D=12 m,厚t=1.0 m,自由端受到y(tǒng)方向的均布載荷p=1 000g(t)。懸臂梁為平面應(yīng)力狀態(tài),計(jì)算參數(shù)取:楊氏模量E=30 MPa,泊松比ν=0.3,質(zhì)量密度ρ=4.0 kg/m3。本例以200個不規(guī)則節(jié)點(diǎn)離散該問題域;采用瑞利阻尼,2個參數(shù)為α和β;應(yīng)用Newmark法求解動力問題方程,取θ=0.5。此時(shí),所采用的方法無條件穩(wěn)定。

圖6 懸臂梁計(jì)算模型

在相同的背景網(wǎng)格下,應(yīng)用2種濃縮形函數(shù)的CSRPIM-T4-Cd和FEM-T3分別計(jì)算該算例,并以有限元軟件用十分精細(xì)的(96×24)四邊形單元計(jì)算的結(jié)果作為參考解,取梁自由端中點(diǎn)A的y方向位移繪制位移時(shí)間歷程曲線,對這3種方法進(jìn)行比較。

取載荷函數(shù)g(t)=cos(ωt)(N/m2),頻率ω=13.5 rad/s,不考慮阻尼的影響(α=0,β=0)。圖7展示的是A點(diǎn)的y方向位移時(shí)間歷程曲線,時(shí)間范圍為2 s,采用小時(shí)間步長(Δt=1 ms)。可以看出,CSRPIM-T4-Cd6的解與參考解十分接近,計(jì)算精度很高。為驗(yàn)證此方法的時(shí)間穩(wěn)定性,在20 s范圍內(nèi)采用較大時(shí)間步長(Δt=5 ms),計(jì)算的位移時(shí)間歷程如圖8所示,可以看出此時(shí)計(jì)算結(jié)果穩(wěn)定,仍然可以得到穩(wěn)定的數(shù)值解。

圖7 余弦載荷下A點(diǎn)在2 s內(nèi)的位移

圖8 CSRPIM-T4-Cd6計(jì)算的余弦載荷下A點(diǎn)在20 s內(nèi)的位移

6 結(jié) 論

本文應(yīng)用濃縮形函數(shù)的單元基光滑點(diǎn)插值法求解二維固體動力學(xué)問題。為驗(yàn)證方法的性質(zhì),進(jìn)行了一系列自由振動和受迫振動的數(shù)值計(jì)算,得到如下結(jié)論:

(1)CSRPIM使用2種形式的濃縮形函數(shù)均能通過標(biāo)準(zhǔn)分片試驗(yàn),且均可以收斂于真實(shí)解;

(2)CSRPIM-T4-Cd6方法能夠準(zhǔn)確地模擬系統(tǒng)剛度,為動力學(xué)問題提供十分精確且穩(wěn)定的解;

(3)CSRPIM-T4-Cd4方法計(jì)算的系統(tǒng)剛度過軟,模態(tài)分析中出現(xiàn)了虛假模態(tài),存在時(shí)間不穩(wěn)定的問題。

(4)虛擬點(diǎn)的加入與濃縮形函數(shù)的使用沒有引入額外的系統(tǒng)自由度,相比傳統(tǒng)的CSRPIM-T4方法不會增加計(jì)算成本。

(5)在可以自動劃分問題域的線性三角形背景網(wǎng)格上,CSRPIM-T4-Cd6方法能夠提供精度很高的計(jì)算結(jié)果,十分適合解決實(shí)際工程問題。