伍洋 邢譽(yù)峰

摘要： 提出了一種基于Hermite插值的微分求積升階譜有限元方法。單元在幾何映射上采用了混合函數(shù)方法，而在形函數(shù)的構(gòu)造上，單元邊界上采用非均勻節(jié)點的Hermite插值基函數(shù)，單元內(nèi)部升階譜形函數(shù)的構(gòu)造則采用雅克比正交多項式的張量積形式。將單元形函數(shù)與高斯-洛巴托積分法結(jié)合起來離散薄板的勢能泛函從而得到相應(yīng)的單元矩陣。提出的薄板單元在單元邊界以及單元內(nèi)部的節(jié)點配置完全自由，因而可以用于不同階次的單元的連接。通過在薄板自由振動中的應(yīng)用計算以及與精確解的對比，結(jié)果表明：提出的微分求積升階譜有限元方法不僅計算精度高，而且收斂速度快，同時在階次較高時仍然具有良好的數(shù)值穩(wěn)定性。

關(guān)鍵詞： 薄板自由振動； 微分求積方法； Hermite插值； 升階譜有限元方法

中圖分類號： O323； O241.82 文獻(xiàn)標(biāo)志碼：A文章編號1004-4523（2018）02-0343-09

DOI：10.16385/j.cnki.issn.1004-4523.2018.02.019

引言

有限元方法自從20世紀(jì)60年代被正式提出以來，目前已經(jīng)廣泛應(yīng)用于工程實際中，并對促進(jìn)經(jīng)濟(jì)發(fā)展，解決工程科學(xué)問題做出了巨大貢獻(xiàn)[1-3]。其中，h-型有限元方法以其技術(shù)簡單、成熟，目前已發(fā)展成為有限元分析的主流方法。通過加密網(wǎng)格，該方法理論上可以得到足夠高精度的計算結(jié)果。然而，在實際計算過程中，h-型有限元方法的收斂速度一般較慢，高精度的計算結(jié)果往往需要非常精細(xì)的網(wǎng)格才能得到，因而對CAD系統(tǒng)的前處理過程要求更加嚴(yán)格。伴隨著h-型方法發(fā)展的還有p-型有限元方法，這種方法在計算過程中可以保持網(wǎng)格劃分不變，通過提高單元的階次便可得到高精度的計算結(jié)果。Babuka 等從理論分析和實際應(yīng)用表明p-型有限元方法具有比h-型方法更快的收斂速度[3-5]。作為一種典型的p-方法，升階譜有限元方法以其獨(dú)特的優(yōu)勢在有限元方法發(fā)展的早期得到了廣泛的研究[6-9]。其最顯著的特點是低階單元的基函數(shù)是高階單元基函數(shù)的子集，這樣在計算過程中中，低階單元的結(jié)構(gòu)矩陣是高階單元結(jié)構(gòu)矩陣的子矩陣，即所謂“嵌套”特性。進(jìn)而在單元階次升高時，前一步的計算結(jié)果仍然可以保留，從而降低重復(fù)分析的成本。基于這些優(yōu)勢，該方法一直被認(rèn)為是未來實現(xiàn)自適應(yīng)分析的主流技術(shù)。然而，由于歷史和該方法自身的技術(shù)困難，升階譜有限元方法的應(yīng)用始終未能得到足夠的推廣[10]。

傳統(tǒng)升階譜有限元方法的基函數(shù)一般不具有插值特性，從而導(dǎo)致節(jié)點自由度不具有明確的物理意義，這樣使得單元的組裝以及邊界條件的施加都變得困難。另外，隨著單元階次升高，高階（幾十階或更高）多項式的計算帶來的數(shù)值穩(wěn)定性問題一直是一個技術(shù)難題[8]。最近，相關(guān)學(xué)者提出了微分求積升階譜有限元方法，并將其應(yīng)用于Mindlin板的自由振動問題中[11]。與傳統(tǒng)升階譜方法不同，該方法在單元邊界上采用了非均勻節(jié)點的插值基函數(shù)，而在單元內(nèi)部則仍采用升階譜面函數(shù)。這樣使得單元組裝和邊界條件施加的問題得到了解決。此外， 該方法采用了正交多項式的迭代公式進(jìn)行計算，有效地克服了高階多項式計算的數(shù)值動蕩特性。

基于微分求積升階譜有限元方法的上述優(yōu)點，本文將把該方法推廣到C1型單元的構(gòu)造上。單元邊界的形函數(shù)的構(gòu)造通過一維高階Hermite插值函數(shù)與三次Hermite插值函數(shù)混合得到，而單元內(nèi)部面函數(shù)的構(gòu)造則通過Jacobi正交多項式的張量積形式得到。其中，單元邊界節(jié)點數(shù)目配置自由，因此可以適用于不同階次的單元拼接以及自適應(yīng)分析中局部自由度的調(diào)整。由于單元邊界插值基函數(shù)的節(jié)點分布將對結(jié)構(gòu)矩陣的數(shù)值特性具有重要的影響，本文提出了一類非均勻節(jié)點，即Gauss-Jacobi點，將有助改善單元的數(shù)值穩(wěn)定性。此外，研究表明這類非均勻節(jié)點具有穩(wěn)定的插值特性，因此可作為微分求積方法（強(qiáng)形式）的非均勻節(jié)點使用。計算表明，本文提出的C1型微分求積升階譜有限元方法具有精度高，收斂速度快等特點，在單元階次較高時，仍然具有較好的數(shù)值穩(wěn)定性。

第2期伍洋，等： 微分求積升階譜有限元方法及其在薄板自由振動中的應(yīng)用振 動 工 程 學(xué) 報第31卷1微分求積升階譜薄板單元〖*2〗1.1一維Hermite插值基函數(shù)對給定定義在區(qū)間[-1，1]上的一維函數(shù)fξ，其Hermite插值近似函數(shù)為fξ≈H0ξf′-1+HN+1ξf′1+

∑Nj=1Hjξfξj=∑N+1j=0Hjξj（1）式中ξj， j=1～N為插值節(jié)點，邊界上的基函數(shù)為H0ξ=1-ξ22L1ξ

H1ξ=c1ξ+c21-ξ2L1ξ

HNξ=c3ξ+c41+ξ2LNξ

HN+1ξ=ξ2-12LNξ（2）內(nèi)部節(jié)點對應(yīng)的基函數(shù)為Hjξ=1-ξ21-ξ2jLjξ， j=2，3，…，N-1（3）式中Lj（ξ），j=1～N，為定義在節(jié)點ξj的拉格朗日插值函數(shù)，而方程中常數(shù)項為：c1=12-L′1ξ1， c2=32-L′1ξ1，

c3=-12-L′NξN， c4=32+L′NξN（4）可以看到內(nèi)部節(jié)點ξj， j=2～N-1對應(yīng)的基函數(shù)Hj （ξ） 滿足在該節(jié)點取值為1，而在其他節(jié)點處取值為0。然而Hj （ξ）在該節(jié)點的取值不一定為最大值，這樣在節(jié)點數(shù)目較多時，其插值精度并不穩(wěn)定，易出現(xiàn)龍格現(xiàn)象。為提高插值精度，令maxHjξ=Hjξj=1， j=2，3，…，N-1（5）這樣得到一組關(guān)于節(jié)點ξj， j=2～N-1的非線性方程組，通過求解上述方程組可以得到一組非均勻分布的節(jié)點。計算表明，該非均勻節(jié)點為權(quán)值為（3， 3）的雅克比正交多項式J（3，3）N-2（ξ） 的0點，因而在本文中記為高斯-雅克比點。文獻(xiàn)[12]給出了多種求解正交多項式0點的方法。圖1給出了當(dāng)N=2和N=10時的Hermite插值基函數(shù)。從圖1（a）中可以看到其明顯的插值特性，從圖1（b）可以看到中間節(jié)點對應(yīng)的基函數(shù)在該節(jié)點取最大值1。

圖1Hermite插值基函數(shù)

Fig.1Hermite interpolation bases1.2幾何映射

在傳統(tǒng)h-型有限元方法中，復(fù)雜結(jié)構(gòu)的幾何形狀可以通過大量單元（通常為直邊單元）來近似表示而不至于引起特別大的誤差。然而，本文所構(gòu)造的單元屬于p-型單元，因此在實際過程中單元的幾何尺寸要比傳統(tǒng)h-型單元尺寸大，無法靠大量的單元來逼近實際的求解域，進(jìn)而要求單元本身能精確地表示幾何形狀。為此，本文將采用混合函數(shù)方法來構(gòu)造精確的曲邊單元，該方法在一般的有限元教材中可以查到[3]。如圖2所示，曲邊單元（右圖）在邊界上的曲線可以通過參考坐標(biāo)系（左圖）表示為

x=xi（ξ），y=yiξ，-1≤ξ≤1，i=1，3

x=xi（η），y=yiη，-1≤η≤1，i=2，4（6）

根據(jù)混合函數(shù)插值，從參考域（母單元）到物理域的幾何映射可以表示為：xξ，η=1-η2x1ξ+1+ξ2x2η+1+η2x3ξ+圖2單元節(jié)點配置與幾何映射

Fig.2Element node collocation and geometric mapping

1-ξ2x4η-1-ξ1-η4X1-

1+ξ1-η4X2-1+ξ1+η4X3-

1-ξ1+η4X4（7）

yξ，η=1-η2y1ξ+1+ξ2y2η+1+η2y3ξ+

1-ξ2y4η-1-ξ1-η4Y1-

1+ξ1-η4Y2-1+ξ1+η4Y3-

1-ξ1+η4Y4（8）式中Xi，Yi，i=1～4為4個頂點的直角坐標(biāo)。此外，在薄板單元中需要計算到二階導(dǎo)數(shù)，因此需要用到如下鏈?zhǔn)椒▌t：x

y=1JJ22-J12

-J21J11ξ

η，

J=J11J12

J21J22=xξyξ

xηyη（9）

2x2

2y2

2xy=1J2·

J222J212-2J12J22

J221J211-2J11J21

-J21J22-J11J12J11J22+J12J21·

2ξ2-2xξ2x-2yξ2y

2η2-2xη2x-2yη2y

2ηξ-2xηξx-2yηξy（10）1.3形函數(shù)

與傳統(tǒng)等參單元類似，本文將在參考坐標(biāo)系下構(gòu)造單元的形函數(shù)，這樣形函數(shù)的參考域與幾何定義的參考域相同，最后在計算過程中通過參數(shù)變換可以得到物理域上的形函數(shù)。如圖2所示，單元在4個頂點各配置4個自由度，分別為w， wξ， wη， wξη， 而在各邊界內(nèi)部每個節(jié)點配置兩個自由度w， wn， 即撓度和法向轉(zhuǎn)角， 各邊上的節(jié)點數(shù)目分別為M， N， P， Q。根據(jù)混合函數(shù)插值，那么可以得到每個自由度對應(yīng)的形函數(shù)，如頂點1SV1w=η+2η-124H（1）1ξ+

ξ+2ξ-124H（1）1η-

η+2η-12ξ+2ξ-1216

SV1wξ=η+2η-124H（1）0ξ+

ξ+1ξ-124H（4）1η-

η+2η-12ξ+1ξ-1216

SV1wη=η+1η-124H（1）1ξ+

ξ+2ξ-124H（4）0η-

η+1η-12ξ+2ξ-1216

SV1wξη=η+1η-124H（1）0ξ+

ξ+1ξ-124H（4）0η-

η+1η-12ξ+1ξ-1216（11）其中上標(biāo)V1代表頂點1，類似地可以得到其他頂點處的形函數(shù)，不再贅述。同理，也可以得到各邊界內(nèi)部節(jié)點對應(yīng)的形函數(shù)，如邊界1上

S（1）w，i=η+2η-124H（1）iξ，i=2～M-1

S（1）wn，i=η+1η-124H（1）iξ，i=2～M-1（12）

式中上標(biāo)（1）代表邊界1，指標(biāo)i代表邊（1）上的第i個節(jié)點。圖3給出了上述形函數(shù)的圖形。

由于薄板單元要滿足C1連續(xù)性，故單元內(nèi)部升階譜形函數(shù)在單元邊界上必須滿足函數(shù)值以及導(dǎo)數(shù)值均為0?？梢酝ㄟ^雅克比正交多項式的張量積形式來進(jìn)行構(gòu)造，即

SFmn=mξnη， m=1～Hξ， n=1～Hη

mξ=ξ2-12J（4，4）m-1ξ（13）

式中J（4，4）m-1為雅克比正交多項式，權(quán)系數(shù)（4， 4）的選擇可以使得上式滿足如下正交性∫1-1∫1-1SFpqSFmndξdη=∫1-1ξ+14ξ-14J（4，4）p-1·

J（4，4）m-1dξ∫1-1η+14η-14J（4，4）q-1J（4，4）n-1dη=

210Γp+42p+7ΓpΓp+8·

210Γm+42m+7ΓmΓm+8δpmδqn（14）形函數(shù)的正交性有利于提高單元結(jié)構(gòu)矩陣的數(shù)值特性。其前兩個形函數(shù)如圖4所示， 可以看到，

圖3 邊界形函數(shù)

Fig.3 Edge shape functions圖4升階譜面函數(shù)

Fig.4Hierarchical face functions面函數(shù)在單元邊界上的函數(shù)值以及導(dǎo)數(shù)值均為0。

1.4數(shù)值離散

薄板的應(yīng)變能與動能系數(shù)可以表示為U=D2S2wx22+2wy22+2ν2wx22wy2+

21-ν2wxy2dxdy

T0=12Sρhw2dxdy（15）式中ν為泊松比，D=Eh3/[12（1-ν2）]為薄板的彎曲剛度，E為彈性模量，ρ為密度，h為板的厚度。將上式離散即可得到單元剛度矩陣以及質(zhì)量矩陣，為此，設(shè)單元的位移試函數(shù)為wξ，η=∑4i=1（SViwwVi+SViwξwViξ+SViwηwViη+

SViwξηwViξη）+∑Hξm=1∑Hηn=1amnSFmn+

∑M-1i=2S（1）w，iw（1）i+S（1）wn，iw（1）n，i+

∑N-1i=2S（2）w，iw（2）i+S（2）wn，iw（2）n，i+

∑P-1i=2S（3）w，iw（3）i+S（3）wn，iw（3）n，i+

∑Q-1i=2S（4）w，iw（4）i+S（4）wn，iw（4）n，i（16）式中wVi，…，wViξη為頂點處的節(jié)點位移；w（j）i，w（j）n，j=1～4為邊內(nèi)的節(jié)點位移；amn為對應(yīng)升階譜面函數(shù)的廣義節(jié)點位移。定義如下節(jié)點位移向量= [wV1，w（1）i（i=2～M-1），

wV2，w（2）i（i=2～N-1），

wV3，w（3）i（i=P-1～2），

wV4，w（4）i（i=Q-1～2），

wV1η，w（1）n，i（i=2～M-1），

wV2η，wV2ξ，w（2）n，i（i=2～N-1），

wV3ξ，wV3η，w（3）n，i（i=P-1～2），

wV4η，wV4ξ，w（4）n，i（i=Q-1～2），

wViξη（i=1～4），a11，a21，…，aHξHη]T通過高斯-洛巴托積分，詳見文獻(xiàn)[11]，式（15）可以離散為如下形式U=12THTBH

T0=12ρhTGTCG（17）其中B=DCυC0

υCC0

0021-υC

H=Gxx

Gyy

Gxy（18）

C=diagC1，C2，…，CNξ，

Ci=CξidiagJi1Cη1，Ji2Cη2，…，JiNηCηNη，

Jij=detJ（ξi，ηj）（19）式中Cξi， Cηj為沿著ξ，η方向的高斯-洛巴托積分系數(shù)；Gxx， Gyy， Gxy為微分矩陣。注意到是參考系下的節(jié)點位移向量，因而不能直接用于單元組裝與邊界條件的施加，為此需要將其轉(zhuǎn)化為x-y坐標(biāo)系中的節(jié)點位移，如圖2（右）所示，為此定義如下節(jié)點位移向量= [wV1，w（1）i（i=2～M-1），

wV2，w（2）i（i=2～N-1），

wV3，w（3）i（i=P-1～2），wV4，

w（4）i（i=Q-1～2），wn，1，…，

wn，M+P+Q+N，wxy，i（i=1～4），

a11，a21，…，aHξHη]T（20）同時可以得到可逆轉(zhuǎn)換關(guān)系=Q（21）將上式代入式（17）可以得到U=12TQ-THTBHQ-1

T0=12ρhTQ-TGTCGQ-1（22）進(jìn)而可以得到剛度矩陣與質(zhì)量矩陣K=Q-THTBHQ-1

M=ρhQ-TGTCGQ-1（23）2算例

本節(jié)將上述單元用于求解薄板的自由振動問題，如圖5所示為一扇形板，其圓弧邊界的半徑分別為a和b， 周向旋轉(zhuǎn)角度為φ。本例中考慮兩種邊界條件，分別為四邊簡支（SSSS） 和兩徑向邊簡支，兩周向邊固支（SCSC）。表1研究了四邊簡支邊界條件下扇形板的前6階固有頻率收斂過程，整個板采用一個單元計算（如圖5所示），節(jié)點數(shù)目在各邊界上相同，記為Ne， 而內(nèi)部兩個方向上的升階譜形函數(shù)的數(shù)目也取為相同值Nh。計算過程中總自由度數(shù)記為TDOFs， 而施加邊界條件之后的自由度數(shù)記為RDOFs。從表中可以看到，隨著邊界節(jié)點數(shù)目的增加，只需要少量的內(nèi)部自由度即可達(dá)到收斂的計算結(jié)果（在表中考慮的有效數(shù)字范圍內(nèi)）。當(dāng)邊界節(jié)點數(shù)目為7，內(nèi)部升階譜面函數(shù)數(shù)目為7×7時，本文計算結(jié)果在考慮的有效數(shù)值范圍內(nèi)即可收斂于精確解，所需要的自由度數(shù)目為105，相比于其他數(shù)值方法，如微分求積有限元方法（DQFEM）、微分求積方法（DQM），所需自由度數(shù)更少。同樣地，對于SCSC邊界條件下的計算結(jié)果也有類似結(jié)論，如表2所示。

圖5扇形板

Fig.5A sectorial plate

為驗證單元組裝時計算特性，考慮圖6所示的正六邊形板，邊長為a/2， 整個板被劃分為3個四邊形單元，各單元邊界上的節(jié)點數(shù)目與內(nèi)部升階譜面函數(shù)關(guān)系為Ne=Nh+2。表3中給出了六邊簡支和六邊固支時的前6階固有頻率的計算結(jié)果和Liew用瑞利-里茲法以及Irie的解析方法結(jié)果。對于六邊簡支情況，由于角點存在彎矩奇異性，因此計算結(jié)果收斂速度較慢，而對于大多數(shù)頻率參數(shù)，計算精度仍然能達(dá)到前3～5位有效數(shù)字收斂，且計算結(jié)果與其他方法相比吻合較好。對于固支邊界情況，速度更快，能達(dá)到前6位有效數(shù)字收斂，且計算結(jié)果均比其他方法給出的頻率要低。

圖6正六邊形板

Fig.6A hexagonal plate表1四邊簡支扇形板的前6階固有頻率收斂過程

Tab.1Convergence of the first 6 natural frequencies of SSSS sectorial plates

a/b=2.0， φ=45°，Ω=ωa2ρh/D

MethodNeNhRDOFs/

TDOFsMode sequence123456Present3757/7368.382152.34190.79282.49286.54390.18872/8868.382152.34190.79282.49286.53390.1720408/42468.382152.34190.79282.49286.53390.174876/9668.380150.99189.61282.42283.76390.17993/11368.380150.99189.61282.41283.76390.1720412/43268.380150.99189.61282.41283.76390.17510116/14068.379150.99189.60278.46283.61387.6611137/16168.379150.99189.60278.46283.61387.6520416/44068.379150.99189.60278.46283.61387.656884/11268.379150.98189.60278.45283.59387.669101/12968.379150.98189.60278.45283.59387.6520420/44868.379150.98189.60278.45283.59387.657660/9268.379150.98189.60278.43283.59387.70773/10568.379150.98189.60278.39283.59387.6220424/45668.379150.98189.60278.39283.59387.62Exact[13]68.379150.98189.60278.39283.59387.62DQFEM[14]N=12—/14468.379150.98189.60278.39——DQM[15]N=21—/44168.380150.96189.60278.35283.54387.62

表2SCSC扇形板的前6階固有頻率收斂過程

Tab.2Convergence of the first 6 natural frequencies of SCSC sectorial plates

a/b=2.0， φ=45°，Ω=ωa2ρh/D

MethodNeNhRDOFs/

TDOFsMode123456Present4985/113107.57178.84269.55306.05346.63477.2610104/132107.57178.84269.55306.05346.63477.2520404/432107.57178.84269.55306.05346.63477.25MethodNeNhRDOFs/

TDOFsMode123456Present5987/121107.57178.83269.54305.86346.56476.3110106/140107.57178.83269.54305.86346.56476.3020406/440107.57178.83269.54305.86346.56476.306989/129107.57178.82269.49305.84346.48476.3110108/148107.57178.82269.49305.84346.48476.3020408/448107.57178.82269.49305.84346.48476.307991/137107.57178.82269.49305.84346.46476.3110110/156107.57178.82269.49305.84346.46476.3020410/456107.57178.82269.49305.84346.46476.30Exact[13]107.57178.82269.49305.84346.46476.30DQFEM[14]N=11-/121107.57178.82269.49305.84——DQM[15]N=21-/441107.57178.79269.50305.80346.40476.24

表3正六邊形板的前6階固有頻率收斂過程

Tab.3Convergence of the first 6 natural frequencies of hexagonal plates

Ω=ωa/π2ρh/D

MethodNe = Nh + 2Mode sequence123456SSSSSSPresent52.9515.7.42297.426313.230813.236015.070772.93167.38887.390513.196113.199115.172892.92117.37407.375013.180813.182615.1970112.91507.36617.366813.172413.173515.2026132.91127.36147.361913.167213.168015.2033152.90867.35847.358813.163813.164415.2027172.90697.35637.356613.161413.161915.2017192.90567.35487.355113.159713.160115.2008212.90467.35377.353913.158513.158715.2001232.90397.35297.353113.157513.157715.1994252.90337.35237.352413.156713.156915.1988272.90297.35177.351913.156113.156315.1983Liew[16]2.92747.39927.415013.213513.240315.2745Irie[17]2.9107.3607.38213.18013.20015.240CCCCCCPresent55.148610.747910.751117.550017.562219.706875.171510.748310.748517.539817.540019.900095.178610.748310.748317.541417.541619.9644115.181210.748310.748317.542117.542219.9879135.182410.748310.748317.542417.542519.9981155.183010.748310.748317.542617.542620.0030175.183310.748310.748317.542717.542720.0057195.183510.748310.748317.542717.542720.0072215.183610.748310.748317.542717.542820.0081235.183610.748310.748317.542817.542820.0087255.183710.748310.748317.542817.542820.0091275.183810.748310.748317.542817.542820.0094Liew[16]5.288610.852510.856717.752817.781620.2331Irie[17]5.21210.85010.88017.73017.75020.130

3結(jié)論

本文提出了一種微分求積升階譜薄板單元，單元形函數(shù)邊界上采用插值型高階Hermite多項式構(gòu)造，其中高斯-雅克比點作為邊界插值節(jié)點有效地提高單元的數(shù)值特性，而內(nèi)部升階譜面函數(shù)采用雅克比正交多項式的張量積來構(gòu)造。該單元在邊界以及內(nèi)部節(jié)點配置自由，單元自由度的配置不隨著網(wǎng)格傳播，可以適用于局部自由度的調(diào)整以及連接不同階次的單元。將本文計算結(jié)果與其他數(shù)值方法以及解析方法的結(jié)果對比表明，該單元不僅計算精度較高，收斂速度快，而且具有良好的數(shù)值穩(wěn)定性。

參考文獻(xiàn)：

[1]Courant R. Variational methods for the solution of problems of equilibrium and vibrations [J]. Bulletin of the American Mathematical Society，1943， 49（1）： 1—23.

[2]Turner M J. Stiffness and deflection analysis of complex structures[J]. Journal of the Aeronautical Sciences， 1956， 23（9）： 805—823.

[3]Zienkiewicz O C， Taylor R L， Taylor R L， et al. Finite Element Method： Its Basis and Fundamentals[M]. Beijing： Beijing World Publishing Corpoaratoin， 2013.

[4]Babuska I， Szabo B A， Katz I N. The p-version of the finite element method[J]. Siam Journal on Numerical Analysis 1981， 18（3）： 515—545.

[5]Babuka I， Suri M.The p- and h-p versions of the finite element method， an overview[J]. Computer Methods in Applied Mechanics and Engineering， 1990， 80（1990）： 5—26.

[6]Zienkiewicz O C， De S R Gago J P， Kelly D W. The hierarchical concept in finite element analysis[J]. Computers & Structures 1983， 16（1）： 53—65.

[7]Peano A. Hierarchies of conforming finite elements for plane elasticity and plate bending[J]. Computers & Mathematics with Applications， 1976， 2（3）： 211—224.

[8]Bardell N S. The application of symbolic computing to the hierarchical finite element method[J]. International Journal for Numerical Methods in Engineering， 1989， 28（5）： 1181—1204.

[9]Bardell N S. Free vibration analysis of a flat plate using the hierarchical finite element method[J]. Journal of Sound & Vibration， 1991， 151（2）： 263—289.

[10]諸德超.升階譜有限元法[M]. 北京：國防工業(yè)出版社， 1993.

[11]Liu C， Liu B， Zhao L， et al. A differential quadrature hierarchical finite element method and its applications to vibration and bending of Mindlin plates with curvilinear domains[J]. International Journal for Numerical Methods in Engineering ， 2016， 109（2）： 174—197.

[12]Shen J， Tang T， Wang L L. Spectral Methods： Algorithms， Analysis and Applications[M]. Berlin： Springer， 2011.

[13]Kim C S， Dickinson S M. On the free， transverse vibration of annular and circular， thin， sectorial plates subject to certain complicating effects[J]. Journal of Sound and Vibration，1989， 134 （3）：407—421.

[14]Xing Y， Liu B. High-accuracy differential quadrature finite element method and its application to free vibrations of thin plate with curvilinear domain[J]. International Journal for Numerical Methods in Engineering， 2009， 80 （13）：1718—1742.

[15]Bert C W， Malik M.The differential quadrature method for irregular domains and application to plate vibration[J]. International Journal of Mechanical Sciences，1996， 38 （6）：589—606.

[16]Liew K M， Lam K Y A. Set of orthogonal plate functions for flexural vibration of regular polygonal plates[J]. Journal of Vibration and Acoustics，1991， 113 （2）：182—186.

[17]Irie T， Yamada G， Umesato K. Free vibration of regular polygonal plates with simply supported edges[J]. The Journal of the Acoustical Society of America， 1981， 69 （5）：1330—1336.

A differential quadrature hierarchical finite element method and its

application to thin plate free vibrationWU Yang， XING Yu-feng（Institute of Solid Mechanics， Beihang University （BUAA）， Beijing 100191， China）

Abstract： A Hermite differential quadrature hierarchical finite element method is proposed and elaborated in this paper. The geometric mapping is based on blending function interpolation， while the shape functions at the element edges are derived from Hermite interpolation bases with non-uniform distributed nodes and the internal hierarchical face functions are obtained through tensor product of one dimensional Jacobi orthogonal polynomials. The shape functions in conjunction with Gauss-Lobatto interpolation are employed in discretizing the potential functional of thin plates to obtain the element matrices. In the present elements， the DOFs （Degree of Freedoms） collocations are free at each side and inside of the element， therefore it is allowed to use elements with different order in assembly. The present elements are used in free vibration analysis of thin plate， and the results are compared with the exact solutions and the numerical results by other methods， which show the high accuracy as well as fast convergence of present elements. Moreover， the stability problems are not suffered even when the element order is very high.

Key words： vibration of thin plates； differential quadrature method； Hermite interpolation； hierarchical finite element method