謝丹 蹇開(kāi)林

摘要： 采用全域插值廣義移動(dòng)最小二乘（IGMLS）的改進(jìn)無(wú)網(wǎng)格伽遼金（EFG）法，對(duì)旋轉(zhuǎn)中心剛體柔性梁系統(tǒng)的剛?cè)狁詈蟿?dòng)力學(xué)特性進(jìn)行研究。利用考慮了柔性梁橫向變形引起的非線性耦合項(xiàng)的一次近似耦合模型，根據(jù)Hamilton 原理和EFG離散方法得到剛?cè)狁詈舷到y(tǒng)的無(wú)網(wǎng)格動(dòng)力學(xué)離散方程。在大范圍運(yùn)動(dòng)未知的情況下，采用數(shù)值方法對(duì)剛?cè)狁詈舷到y(tǒng)進(jìn)行動(dòng)力響應(yīng)的仿真計(jì)算，并對(duì)EFG法的主要影響因素進(jìn)行了討論分析。通過(guò)與有限元法的數(shù)值計(jì)算結(jié)果對(duì)比，驗(yàn)證EFG法用于剛?cè)狁詈舷到y(tǒng)動(dòng)力學(xué)研究的有效性及可行性，并為無(wú)網(wǎng)格法用于更復(fù)雜的柔性多體動(dòng)力學(xué)的研究提供了理論依據(jù)。關(guān)鍵詞： 多體動(dòng)力學(xué)； 剛?cè)狁詈希?旋轉(zhuǎn)梁； 無(wú)網(wǎng)格伽遼金（EFG）法； 廣義移動(dòng)最小二乘（GMLS）

中圖分類號(hào)：O313.7文獻(xiàn)標(biāo)志碼： A文章編號(hào)： 10044523（2017）04052708

DOI：10.16385/j.cnki.issn.10044523.2017.04.001

引言

作為柔性多體動(dòng)力學(xué)的一種典型模型，旋轉(zhuǎn)梁的動(dòng)力學(xué)特性已被廣泛而深入的研究，包括附著在運(yùn)動(dòng)基上的結(jié)構(gòu)動(dòng)力學(xué)問(wèn)題（或動(dòng)力剛化問(wèn)題）以及大范圍運(yùn)動(dòng)情況未知的剛?cè)狁詈蟿?dòng)力學(xué)問(wèn)題[17]。在柔性多體動(dòng)力學(xué)的研究中，柔性體變形場(chǎng)的離散方法是其中的一個(gè)關(guān)鍵問(wèn)題，根據(jù)目前的國(guó)內(nèi)外研究，對(duì)柔性體變形場(chǎng)的離散主要采用的是假設(shè)模態(tài)法和有限元法[4，5]。其中，假設(shè)模態(tài)法由于依賴已知結(jié)構(gòu)的振動(dòng)振型而使其應(yīng)用往往限于簡(jiǎn)單結(jié)構(gòu)，有限元法雖具有很強(qiáng)的通用性，其表現(xiàn)出的固有缺陷則源于對(duì)單元或網(wǎng)格的依賴。

近年來(lái)，無(wú)網(wǎng)格法[8]以其獨(dú)特的優(yōu)勢(shì)得到了國(guó)內(nèi)外各領(lǐng)域?qū)W者的研究。它通過(guò)一組無(wú)需事先定義、散布在問(wèn)題域及邊界上的節(jié)點(diǎn)來(lái)表示問(wèn)題域及邊界，然后用場(chǎng)節(jié)點(diǎn)來(lái)構(gòu)造近似的總體函數(shù)，可以徹底或部分的消除網(wǎng)格。有限元法采用預(yù)定義的單元構(gòu)造形函數(shù)，所有單元的形函數(shù)是相同的；而無(wú)網(wǎng)格法的近似函數(shù)采用的是定義在局部支持域中的緊支函數(shù)，任意計(jì)算點(diǎn)的形函數(shù)僅通過(guò)包含在其支持域內(nèi)的場(chǎng)節(jié)點(diǎn)來(lái)構(gòu)造，計(jì)算點(diǎn)處的形函數(shù)可隨計(jì)算點(diǎn)的變化而變化。因此，形函數(shù)的構(gòu)造差異是無(wú)網(wǎng)格法與有限元法的本質(zhì)區(qū)別。目前主要的無(wú)網(wǎng)格法有：無(wú)網(wǎng)格伽遼金（EFG）法[9]、徑向點(diǎn)插值法（RPIM）[10]及再生核粒子法（RKPM）[11]等。無(wú)網(wǎng)格法的研究雖已被用于各種領(lǐng)域[1216]，但無(wú)網(wǎng)格法在多體動(dòng)力學(xué)中的研究鮮有報(bào)道。文獻(xiàn)[17]采用點(diǎn)插值法對(duì)旋轉(zhuǎn)懸臂梁的彎曲振動(dòng)及非慣性系下的結(jié)構(gòu)動(dòng)力學(xué)特性進(jìn)行了分析，忽略了大范圍運(yùn)動(dòng)對(duì)變形運(yùn)動(dòng)的影響。

眾多無(wú)網(wǎng)格法中，EFG法以其高精度及簡(jiǎn)明的變換關(guān)系得到了廣泛的研究與應(yīng)用[1516]，但傳統(tǒng)EFG法由于基于移動(dòng)最小二乘近似（MLS）[18]，位移及導(dǎo)數(shù)邊界的施加比較困難[1920]。本文采用全域插值的廣義移動(dòng)最小二乘近似（IGMLS）對(duì)柔性梁的變形場(chǎng)進(jìn)行離散，在保持EFG法高精度的同時(shí)，使其位移及導(dǎo)數(shù)邊界條件的施加像有限元法一樣直接簡(jiǎn)單。在此基礎(chǔ)上，采用考慮了非線性耦合項(xiàng)的一次近似模型[5]，對(duì)大范圍運(yùn)動(dòng)未知情況下的旋轉(zhuǎn)中心剛體柔性梁系統(tǒng)的剛?cè)狁詈蟿?dòng)力學(xué)特性進(jìn)行了數(shù)值研究，并通過(guò)與有限元法對(duì)比，論證了EFG法用于剛?cè)狁詈蟿?dòng)力學(xué)研究的可行性。

1剛?cè)狁詈舷到y(tǒng)運(yùn)動(dòng)學(xué)及變形描述

假定系統(tǒng)中的柔性梁采用的是Euler梁模型，材料均勻且各向同性，梁上的各處的橫截面積均相等，不計(jì)重力影響。

圖1為中心剛體懸臂梁系統(tǒng)，中心剛體繞Z0軸轉(zhuǎn)動(dòng)，研究系統(tǒng)的剛體運(yùn)動(dòng)和懸臂梁在O0X0Y0平面內(nèi)的運(yùn)動(dòng)與變形。其中： O0X0Y0Z0為固定坐標(biāo)系（右手定則），其原點(diǎn)與中心剛體的回轉(zhuǎn)中心重合，OXY為固結(jié)在梁端點(diǎn)的浮動(dòng)坐標(biāo)系。梁的中性軸線上任意一點(diǎn)P0發(fā)生變形后的位置為P，rA為浮動(dòng)坐標(biāo)系原點(diǎn)O在固定坐標(biāo)系下的位置矢量，h為P0在浮動(dòng)坐標(biāo)系下的位置矢量，u為P相對(duì)P0的變形位移矢量，r為P在固定坐標(biāo)系下的位置矢量。相關(guān)幾何尺寸及物理參數(shù)表示如下：a為中心剛體的半徑，JH為中心剛體關(guān)于Z0軸的轉(zhuǎn)動(dòng)慣量，L為柔性梁未變形前的長(zhǎng)度，ρ為梁的體積密度，A為梁的橫截面積，E為材料的彈性模量，I為梁的橫截面關(guān)于中性軸（Z軸）的慣性矩，τ為作用在中心剛體上的驅(qū)動(dòng)力矩（力偶矩矢量沿Z0軸），θ為中心剛體作大范圍旋轉(zhuǎn)運(yùn)動(dòng)的角位移。

圖1旋轉(zhuǎn)中心剛體柔性梁系統(tǒng)物理模型

Fig.1The physical model of rotatinghub beam system

第4期謝丹，等：改進(jìn)EFG法用于旋轉(zhuǎn)梁的剛?cè)狁詈蟿?dòng)力學(xué)研究振 動(dòng) 工 程 學(xué) 報(bào)第30卷在慣性坐標(biāo)系下，梁中性軸線上任意一點(diǎn)P的位置矢量r對(duì)應(yīng)的坐標(biāo)向量為r=rA+Θ（h+u）（1）其中：h=x，0T （2）

Θ=cosθ-sinθ

sinθcosθ （3）

rA=acosθ

asinθ=aΘe0 （4）

u=u1

u2=ω1+ωc

ω2（5）式中e0=1，0T，式（5）為一次近似耦合模型的變形位移表達(dá)式，ω1與ω2分別為梁中軸線上一點(diǎn)沿X軸的位移和沿Y軸的位移（即橫向變形），ωc=-12∫x0（ω2x）2dx為梁橫向變形引起的軸向縮短量，即耦合變形量[5]，若不考慮該耦合變形量，則為零次近似模型。

對(duì)式（1）求導(dǎo)，得到點(diǎn)P的速度向量（絕對(duì)速度在固定坐標(biāo)系下的投影）=A+Θ（h+u）+Θ （6）對(duì)式（6）求導(dǎo)，得到點(diǎn)P的加速度矢量（絕對(duì)加速度在固定坐標(biāo)系下的投影）=A-2Θ（h+u）+2Θ+Θ+

Θ（h+u）（7）式中=0-1

10（8）系統(tǒng)動(dòng)能的變分δT=-JHδθ-∫L0ρAδrTdx （9）考慮非線性應(yīng)變位移關(guān)系式εxx=u1x-y2u2x2+12（u2x）2，得到勢(shì)能的變分δΠ=∫L0EAω1xδ（ω1x）dx+∫L0EI2ω2x2δ2ω2x2dx（10）作用在系統(tǒng)上的外力的虛功δWF=τδθ（11）2系統(tǒng)剛?cè)狁詈蟿?dòng)力學(xué)離散方程2.1全域插值廣義移動(dòng)最小二乘（IGMLS）近似假設(shè)歐拉梁的求解域0≤x≤L被用N個(gè)場(chǎng)節(jié)點(diǎn)xI（I=1， 2，…，N）離散，根據(jù)EFG方法，撓度變量v（x）在全求解域內(nèi)的近似表達(dá)式為vh（x）=∑mj=1pj（x）aj=pT（x）a （12）式中a=a1，a2，…，amT為待定參數(shù)向量， p（x）=1，x，x2T為m=3時(shí)的基函數(shù)向量。

則撓度函數(shù)及一階導(dǎo)數(shù)（轉(zhuǎn)角）的近似值在各節(jié)點(diǎn)處誤差的加權(quán)平方和為J（a）=∑NI=1ωI（x）vh（xI）-vI2+

θh（xI）-θI2=

∑NI=1ωI（x）∑mj=1pj（x）aj-vI2+

∑mj=1dpj（x）dxaj-θI2（13）式（13）中的權(quán)系數(shù)ωI（x）=ω（x-xI），以下為本文采用的三次樣條權(quán)函數(shù)的表達(dá)式ωI（r）=23-4r2+4r3，r≤0.5

43-4r+4r2+4r33，0.5

0，r>1（14）式中r=dIdmI，dI=‖x-xI‖是計(jì)算點(diǎn)x到節(jié)點(diǎn)xI的距離，節(jié)點(diǎn)影響域尺寸dmI=scale×cI，cI通常取節(jié)點(diǎn)xI與其最近的節(jié)點(diǎn)之間的距離，結(jié)構(gòu)動(dòng)力學(xué)分析中，影響域因子scale通常在2～4取值[19]。

根據(jù)GMLS近似，式（13）中的待定參數(shù)要使誤差平方和最小，計(jì)算整理得A（x）a=B（x）（15）式中=1，2，…，N，1，2，…，NT（16）

A（x）=∑NI=1ωI（x）p（xI）pT（xI）+

dp（xI）dxdpT（xI）dx（17）

B（x）=B（0）（x）B（1）（x）

B（0）（x）=ω1（x）p（x1）…ωN（x）p（xN）

B（1）（x）=ω1（x）dp（x1）dx…ωN（x）dp（xN）dx （18） 將式（15）代入（12）得到梁的GMLS形函數(shù)，即Ω（1×2N）=ΦΨ=pT（x）A（x）-1B（x）（19）因此，基于節(jié)點(diǎn)位移參數(shù)的撓度及轉(zhuǎn)角近似表達(dá)式為vh（x）=∑NI=1IvI+∑NI=1ΛI(xiàn)θI=Ω

θh（x）=∑NI=1I，xvI+∑NI=1ΛI(xiàn)，xθI=Ω，x （20）真實(shí)節(jié)點(diǎn)位移矢量d與節(jié)點(diǎn)位移參數(shù)矢量的關(guān)系d=Λ2（21）式中Λ2為轉(zhuǎn)換矩陣，Λ（2N×2N） =

1 （x1 ） … N （x1 ） ψ1 （x1 ） … ψN （x1 ）

1 （xN ） … N （xN ） ψ1 （xN ） … ψN （xN ）

1，x （x1 ） … N，x （x1 ） ψ1，x （x1 ） … ψN，x （x1 ）

1，x （xN ） … N，x （xN ） ψ1，x （xN ） … ψN，x （xN ）（22）將式（21）代入式（20）得vh（x） = Ω =ΩΛ-12 d =N2d

θh（x） =Ω，x =Ω，x Λ-12 d= N2，x d（23）則梁的IGMLS形函數(shù)及其一階導(dǎo)數(shù)表達(dá)式如下N2 =ΩΛ-12

N2，x =Ω，x Λ-12 （24）圖2給出了梁離散后，中間節(jié)點(diǎn)的撓度與轉(zhuǎn)角變量對(duì)應(yīng)的標(biāo)準(zhǔn)GMLS形函數(shù)及IGMLS形函數(shù)曲線，由圖可以明顯看出IGMLS的全域插值特性。圖2梁中間節(jié)點(diǎn)的GMLS形函數(shù)及IGMLS形函數(shù)

Fig.2Shape functions of GMLS and IGMLS at the middle point of the beam

2.2剛?cè)狁詈蟿?dòng)力學(xué)離散方程

假設(shè)N1（x）∈R1×N表示軸向位移的插值性移動(dòng)最小二乘形函數(shù)（IMLS），N2（x）∈R1×2N表示上節(jié)推導(dǎo)的橫向位移的IGMLS形函數(shù)，q1（t）∈RN表示節(jié)點(diǎn)縱坐標(biāo)列向量，q2（t）∈R2N表示由撓度與轉(zhuǎn)角組成的節(jié)點(diǎn)坐標(biāo)列向量，表示如下q1（t）=u1（t），u2（t），…，uN（t）T

q2（t）=v1（t），…，vN（t），θ1（t），…，θN（t）T（25）梁的軸向、橫向變形量及耦合變形量的離散表達(dá)式分別為：ω1（x，t）=N1（x）q1（t）（26）

ω2（x，t）=N2（x）q2（t）（27）

ωc （x，t） = -12qT2 S（x）q2 （28）式（28）中，S（x）為耦合形函數(shù)S（x） = ∫x0（N2 ′）TN2 ′dx （29）式中N2 ′表示對(duì)坐標(biāo)x求偏導(dǎo)。

根據(jù)Hamilton變分原理，將相關(guān)離散變量代入，省略相關(guān)高階項(xiàng)，考慮到δθ，δqT1 ，δqT2 的任意性，得到系統(tǒng)的剛?cè)狁詈蟿?dòng)力學(xué)無(wú)網(wǎng)格離散方程MθθMθq1Mθq2

Mq1θMq1q10

Mq2θ0Mq2q2

1

2+

2000

00Gq1q2

0Gq2q10

1

2+

000

0Kq1q10

00Kq2q2θ

q1

q2=

Qθ

Qq1

0+τ

0

0（30）式中：Mθθ = JH + Jb +qT1 M1 q1 +

qT2 M2 q2 + 2U11 q1 -qT2 D1 q2 （31）

Mq1θ=（Mθq1）T=-Rq2（32）

Mq2θ=（Mθq2）T=（U12）T+RTq1（33）

Mq1q1=M1（34）

Mq2q2=M2 （35）

Gq1q2=-（Gq2q1）T=-R （36）

Kq1q1=K1-2M1（37）

Kq2q2=K2+2（D1-M2）（38）

Qθ = -2（qT1 M1 1 + qT2 M2 2 +

U11 1 -qT2 D1 2 ） （39）

Qq1 = 2UT11 （40）式（31）中Jb=∫L0ρA（a+x）2dx為柔性梁的轉(zhuǎn)動(dòng)慣量。

式（31）～（40）中常數(shù)矩陣采用背景網(wǎng)格積分，以下給出部分矩陣的EFG法數(shù)值計(jì)算表達(dá)式：M1 = ∫L0ρANT1 N1 dx =

∑nck = 1∑ngi = 1ρAN1 （xQi ）TN1 （xQi ）iJik （41）

D1=∫L0ρA（a+x）S（x）dx=

∑nck=1∑ngi=1ρA（a+xQi）S（xQi）iJik（42）式中nc為用于積分的背景網(wǎng)格數(shù)目，ng 是每個(gè)背景網(wǎng)格內(nèi)的Gauss積分點(diǎn)數(shù)，xQi為計(jì)算點(diǎn)坐標(biāo)，i為第i個(gè)Gauss積分點(diǎn)的加權(quán)系數(shù)，Jik為第k個(gè)背景網(wǎng)格在計(jì)算點(diǎn)xQi處積分的雅克比矩陣。

3數(shù)值算例

對(duì)于剛?cè)狁詈蟿?dòng)力系統(tǒng)，通常外部作用力或力矩規(guī)律是已知的，而大范圍運(yùn)動(dòng)量和梁的振動(dòng)變形量均為求解量。不考慮阻尼影響，柔性梁的物理參數(shù)：長(zhǎng)度L=1.8 m，橫截面積A=2.5×10-4 m2，材料密度ρ=2.7667×103 kg/m3，彈性模量E=6.8952×1010 Pa，橫截面的慣性矩I=1.3021×10-10m4。取中心剛體半徑為a=0.05 m，轉(zhuǎn)動(dòng)慣量JH=0.3 kg·m2，作用在剛體上的外力矩為τ（t）=τ0sin（2πTt），0≤t≤T

0，t> T （43）式中T=2 s，按τ0=1 N·m和τ0=7 N·m兩種情況進(jìn)行數(shù)值仿真計(jì)算。

數(shù)值分析中，采用Newmark法對(duì)結(jié)構(gòu)的非線性動(dòng)力響應(yīng)進(jìn)行計(jì)算，時(shí)間步長(zhǎng)取0.05 s，Newmark參數(shù)β=3/4，γ=1/2。有限元法采用9個(gè)單元離散，無(wú)網(wǎng)格法采用10個(gè)離散節(jié)點(diǎn)，且影響域因子取2.5，選擇四點(diǎn)高斯積分。

3.1剛?cè)狁詈蟿?dòng)力響應(yīng)分析

用FEM1表示采用一次近似模型的有限元法計(jì)算結(jié)果，EFG1表示采用一次近似模型的無(wú)網(wǎng)格法計(jì)算結(jié)果，EFG0則為采用零次近似模型的無(wú)網(wǎng)格法計(jì)算結(jié)果。

圖3～6為τ0=1 N·m時(shí)的數(shù)值計(jì)算結(jié)果， 圖3，4分別為剛體的角位移響應(yīng)與角速度響應(yīng)，圖5，6分別為梁末端的橫向、縱向位移響應(yīng)。由圖可看出，不論大范圍旋轉(zhuǎn)運(yùn)動(dòng)還是柔性體的變形運(yùn)動(dòng)，一次近似耦合模型下的無(wú)網(wǎng)格法與有限元法計(jì)算結(jié)果高

圖3中心剛體轉(zhuǎn)動(dòng)角位移（τ0=1 N·m）

Fig.3The angular displacement of the hub when τ0=1 N·m圖4中心剛體的轉(zhuǎn)速（τ0=1 N·m）

Fig.4The angular velocity of the hub when τ0=1 N ·m圖5柔性梁末端的橫向位移（τ0=1 N·m）

Fig.5The transverse displacement of the tip of the beam when τ0=1 N·m圖6柔性梁末端的縱向位移（τ0=1 N·m）

Fig.6The axial displacement of the tip of the beam when τ0=1 N·m度逼近。當(dāng)外力矩撤銷后，中心剛體大約在20.77°附近擺動(dòng)，柔性梁末端的橫向振動(dòng)頻率約為1.31 Hz，柔性梁的殘余振動(dòng)與中心剛體的轉(zhuǎn)動(dòng)相互激勵(lì)，由于未考慮阻尼影響，相互影響的兩種運(yùn)動(dòng)將一直保持下去，并且在變形運(yùn)動(dòng)中，梁末端的縱向位移遠(yuǎn)遠(yuǎn)小于橫向位移。此外，不考慮耦合效應(yīng)的零次模型在此情況下的計(jì)算結(jié)果也比較接近。

圖7～10為τ0=7 N·m時(shí)的數(shù)值計(jì)算結(jié)果。與τ0=1 N·m時(shí)的計(jì)算結(jié)果相比，由于驅(qū)動(dòng)力矩的增大，剛體的轉(zhuǎn)速明顯提升，柔性梁橫、縱向位移的幅值也明顯增大。在τ0=7 N·m下，基于一次模型的無(wú)網(wǎng)格法與有限元法的計(jì)算結(jié)果依然高度逼圖7中心剛體轉(zhuǎn)動(dòng)角位移（τ0=7 N·m）

Fig.7The angular displacement of the hub when τ0=7 N·m圖8中心剛體的轉(zhuǎn)速（τ0=7 N·m）

Fig.8The angular velocity of the hub when τ0=7 N·m圖9柔性梁末端的橫向位移（τ0=7 N·m）

Fig.9The transverse displacement of the tip of the beam when τ0=7 N·m

圖10柔性梁末端的縱向位移（τ0=7 N·m）

Fig.10The axial displacement of the tip of the beam when τ0=7 N·m近，中心剛體大約在147.1°附近擺動(dòng)，柔性梁末端的橫向振動(dòng)頻率約為1.31 Hz，大范圍運(yùn)動(dòng)與變形運(yùn)動(dòng)的耦合效應(yīng)明顯，而零次模型則表現(xiàn)出很大的差異。進(jìn)一步研究發(fā)現(xiàn)，隨著驅(qū)動(dòng)力矩增大，一次模型下梁末端的位移響應(yīng)可以始終保持穩(wěn)定，而零次模型在τ0=10 N·m時(shí)，數(shù)值結(jié)果發(fā)散。

此外，由圖7～10發(fā)現(xiàn)，零次模型下響應(yīng)的結(jié)果不僅與一次模型有較大差異，而且其幅值相對(duì)一次模型的幅值偏小，這是由于該數(shù)值算例中所取的中心剛體的轉(zhuǎn)動(dòng)慣量（JH=0.3 kg·m2）遠(yuǎn)小于柔性梁的轉(zhuǎn)動(dòng)慣量（Jb≈1.46 kg·m2）的緣故。進(jìn)一步研究發(fā)現(xiàn)，在零次模型發(fā)散之前，隨著驅(qū)動(dòng)力矩增大，當(dāng)中心剛體轉(zhuǎn)動(dòng)慣量遠(yuǎn)小于柔性梁的轉(zhuǎn)動(dòng)慣量（JHJb） 時(shí)，零次模型下系統(tǒng)響應(yīng)的幅值小于一次模型下的幅值（如圖11所示）；當(dāng)中心剛體轉(zhuǎn)動(dòng)慣量與柔性梁的轉(zhuǎn)動(dòng)慣量相當(dāng)時(shí)（JH≈Jb），零次模型下的計(jì)算結(jié)果則偏大（如圖12所示）；當(dāng)中心剛體轉(zhuǎn)動(dòng)慣量遠(yuǎn)大于柔性梁的轉(zhuǎn)動(dòng)慣量（JHJb），零次模型與一次模型的響應(yīng)幅值相近（如圖13所示）。

圖11不同驅(qū)動(dòng)力偶矩下柔性梁末端的橫向位移幅值（JHJb）

Fig.11The amplitude of transverse displacements of the beam′s tip when JHJb圖12不同驅(qū)動(dòng)力偶矩下柔性梁末端的橫向位移幅值（JH≈Jb）

Fig.12The amplitude of transverse displacements of the beam′s tip when JH≈Jb

圖13不同驅(qū)動(dòng)力偶矩下柔性梁末端的橫向位移幅值（JHJb）

Fig.13The amplitude of transverse displacements of the beam′s tip when JHJb3.2影響因素分析

EFG法中的影響域因子及Gauss積分點(diǎn)數(shù)目對(duì)計(jì)算結(jié)果的影響分別如表1和2所示。表中，把FEM的計(jì)算耗時(shí)記作1，EFG法與FEM法的計(jì)算耗時(shí)的比值為計(jì)算相對(duì)耗時(shí)，為說(shuō)明計(jì)算的效果，選擇2 s時(shí)刻梁末端的橫向位移的計(jì)算結(jié)果列于表中，位移相對(duì)誤差則通過(guò)EFG-FEMFEM%計(jì)算得到。

由表1可以看出，隨著影響域因子的增大，EFG法的計(jì)算耗時(shí)逐漸增多，但在所取的影響域因子范圍內(nèi)，EFG法的耗時(shí)小于有限元法。從2 s時(shí)刻梁末端的橫向位移的計(jì)算誤差來(lái)看，影響域因子對(duì)計(jì)算精度的影響不是很明顯，也沒(méi)有表現(xiàn)出特定的規(guī)律，說(shuō)明EFG法在所取影響域因子范圍內(nèi)計(jì)算比較穩(wěn)定，且在達(dá)到與有限元法相似的精度下，其計(jì)算效率高于有限元法。表1不同影響域因子下EFG法的計(jì)算相對(duì)耗時(shí)及橫向位移誤差

Tab.1The relative computational cost and computational error of EFG method based on different influence factors

離散

方式影響域

因子計(jì)算相對(duì)

耗時(shí)2 s時(shí)刻的末端

橫向位移/m

（相對(duì)誤差/%）FEM——1-0.0195（——）1.50.250-0.0196（0.513）2.00.296-0.0195（0.0）EFG2.50.312-0.0194（0.513）3.00.343-0.0186（4.615）4.00.530-0.0195（0.0）

表2 不同Gauss積分點(diǎn)數(shù)下EFG法的計(jì)算相對(duì)耗時(shí)及橫向位移誤差

Tab.2The relative computational cost and computational error of EFG method based on different Gauss integral points

離散

方式Gauss積分

點(diǎn)數(shù)目計(jì)算相對(duì)

耗時(shí)2 s時(shí)刻的末端

橫向位移/m

（相對(duì)誤差/%）FEM——1-0.0195（——）20.2340.054（376.92）EFG30.296-0.0195（0.0）40.310-0.0195（0.0）60.359-0.0172（11.79）70.390-0.0195（0.0）

由于近似函數(shù)構(gòu)造方法的原因，無(wú)網(wǎng)格法的被積函數(shù)通常是有理式，一般情況下無(wú)法進(jìn)行精確積分，采用Gauss積分的積分點(diǎn)數(shù)對(duì)計(jì)算結(jié)果有明顯影響。如表2所示，隨著積分點(diǎn)數(shù)的增多，計(jì)算耗時(shí)逐漸增加，不恰當(dāng)?shù)姆e分點(diǎn)的選取，還會(huì)造成很大的誤差，如積分點(diǎn)取2和6時(shí)。通過(guò)本文的數(shù)值試驗(yàn)，通常Gauss積分點(diǎn)取4即可滿足計(jì)算效率和精度的要求。

4結(jié)論

1）傳統(tǒng)EFG法雖具有高的計(jì)算精度，但邊界條件施加困難；改進(jìn)后的EFG法通過(guò)采用全域插值的廣義移動(dòng)最小二乘（IGMLS）近似，在保持高計(jì)算精度的同時(shí)，使中心剛體旋轉(zhuǎn)梁系統(tǒng)的位移及導(dǎo)數(shù)邊界的施加像有限元法一樣方便處理。

2）零次模型在作用力矩增大到一定值時(shí)出現(xiàn)發(fā)散，而一次模型保持穩(wěn)定。EFG法證明了零次模型在剛?cè)狁詈戏治鲋械木窒扌砸约耙淮谓颇Ｐ偷暮侠硇浴?/p>

3）EFG法中的影響域因子在1.0

參考文獻(xiàn)：

[1]Kane T R， Ryan R R， Banerjee A K， Dynamics of a cantilever beam attached to a moving base [J]. Journal of Guidance， Control and Dynamics， 1987， 10（2）：139—151.

[2]Yoo H H， Ryan R R， Scott R A. Dynamics of flexible beams undergoing large overall motions [J]. Journal of Sound and Vibration， 1995， 181： 261—278.

[3]劉錦陽(yáng)，洪嘉振. 剛?cè)狁詈蟿?dòng)力學(xué)系統(tǒng)的建模理論研究[J].力學(xué)學(xué)報(bào)， 2002， 34（3）：409—415.

Liu Jinyang， Hong Jiazhen. Study on dynamic modeling theory of rigidflexible coupling systems [J]. Acta Mechanica Sinica 2002， 34（3）： 409—415.

[4]蔡國(guó)平， 洪嘉振. 旋轉(zhuǎn)運(yùn)動(dòng)柔性梁的假設(shè)模態(tài)方法研究[J].力學(xué)學(xué)報(bào)， 2005， 37（1）：48—56.

Cai Guoping， Hong Jiazhen. Assumed mode method of a rotating flexible beam [J]. Acta Mechanica Sinica 2005， 37（1）： 48—56.

[5]Yang H， Hong J Z， Yu Z Y， Dynamics modeling of a flexible hubbeam system with a tip mass[J]. Journal of Sound and Vibration， 2003， 266：759—774.

[6]章定國(guó)， 余紀(jì)邦. 做大范圍運(yùn)動(dòng)的柔性梁的動(dòng)力學(xué)分析[J].振動(dòng)工程學(xué)報(bào)， 2006， 19（4）：475—480.

Zhang Dingguo， Yu Jibang. Dynamical analysis of a flexible cantilever beam with large overall motions [J]. Journal of Vibration Engineering 2006， 19（4）： 475—480.

[7]吳勝寶， 章定國(guó). 大范圍運(yùn)動(dòng)剛體柔性梁剛?cè)狁詈蟿?dòng)力學(xué)分析[J].振動(dòng)工程學(xué)報(bào)， 2011， 24（1）：1—7.

Wu Shengbao， Zhang Dingguo. Rigidflexible coupling dynamic analysis of hubflexible beam with large overall motion [J]. Journal of Vibration Engineering 2011， 24（1）： 1—7.

[8]Belytschko T， Kronganz Y， Organ D， et al. Meshless methods： an overview and recent developments [J]. Computer Methods in Applied Mechanics and Engineering， 1996，139（1） ： 3—47.

[9]Belytschko T， Lu Y Y， Gu L. Elementfree Galerkin methods [J]. International Journal for Numerical Methods in Engineering， 1994， 37： 229—256.

[10]Gu Y T， Liu G R. A local point interpolation method for static and dynamic analysis of thin beams [J]. Computer Methods in Applied Mechanics and Engineering， 2001，190： 5515—5528.

[11]Liu W K， Jun S， Zhang Y F， et al. Reproducing kernel particle methods for structural dynamics[J]. International Journal for Numerical Methods in Engineering， 1995， 38： 1655—1679.

[12]Liu L， Chua L P， Ghista D N. Elementfree Galerkin method for static and dynamic analysis of spatial shell structures [J]. Journal of Sound Vibration， 2006， 295： 388—406.

[13]Singh I V， Tanaka M， Endo M. Thermal analysis of CNTbased nanocomposites by element free Galerkin method [J]. Computational Mechanics， 2007， 39： 719—728.

[14]Mirzaei D， Schaback R. Solving heat conduction problems by the Direct Meshless Local PetrovGalerkin （DMLPG） method [J]. Numerical Algorithms， 2014， 65： 275—291.

[15]Hajiazizi M， Bastan P. The elastoplastic analysis of a tunnel using the EFG method： A comparison of the EFGM with FEM and FDM [J]. Applied Mathematics and Computation， 2014， 234： 82—113.

[16]Zhang L W， Deng Y J， Liew K M. An improved element—free Galerkin method for numerical modeling of the biological population problems[J]. Engineering Analysis with Boundary Elements， 2014， 40：181—188.

[17]杜超凡， 章定國(guó)， 洪嘉振. 徑向基點(diǎn)插值法在旋轉(zhuǎn)柔性梁動(dòng)力學(xué)中的應(yīng)用[J]. 力學(xué)學(xué)報(bào)， 2015， 47（2）：279—288.

Du Chaofan， Zhang Dingguo， Hong Jiazhen. A meshfree method based on radial point interpolation method for the dynamic analysis of rotating flexible beams [J]. Chinese Journal of Theoretical and Applied Methanics 2015， 47（2）： 279—288.

[18]張雄， 劉巖. 無(wú)網(wǎng)格法[M]. 北京： 清華大學(xué)出版社， 2004.

[19]Cheng Y M， Bai F N， Peng M J. A novel interpolating elementfree Galerkin （IEFG） method for twodimensional elastoplasticity [J]. Applied Mathematical Modelling， 2014， 38：5187—5197

[20]Garijo D， Valenciaó F， GómezEscalonilla F J. Global interpolating MLS shape functions for structural problems with discrete nodal essencial boundary conditions [J]. Acta Mechanica， 2015， 226： 2255—2276.

An improved EFG approach for rigidflexible coupling dynamics

of the rotating hubbeam system

XIE Dan1， JIAN Kailin1，2

（1. College of Aerospace Engineering， Chongqing University， Chongqing 400044， China；

2. State Key Laboratory of Coal Mine Disaster Dynamics and Control， Chongqing 400044， China）

Abstract： Based on the global interpolating generalized moving least squares （IGMLS）， the improved elementfree Galerkin （EFG） method was applied to rigidflexible coupling dynamics of rotating hubbeam systems. With the firstorder approximate model including the nonlinear coupling term induced by the transverse deformation of the flexible beam， the elementfree dynamic equations were derived by the Hamilton principle and EFG method. Research on the rigidflexible coupling dynamics where the overall motion was unknown was conducted by numerical methods； meanwhile， the major influence factors of the EFG method were discussed. Numerical comparisons between the EFG method and the FEM verified the effectiveness and feasibility of the EFG method applied for rigidflexible coupling dynamics and this research provided a direction of elementfree methods for more complicated flexible multibody dynamics.Key words： multibody dynamics； rigidflexible coupling； rotating beam；elementfree Galerkin （EFG） method；generalized moving least squares （GMLS）作者簡(jiǎn)介：謝丹（1987—），女，博士研究生。電話：158238266528；Email： xiedanv@163.com

通訊作者：蹇開(kāi)林（1965—），男，博士，教授。電話：13220275685；Email：cqjian@cqu.edu.cn