微分求積升階譜有限元方法及其在薄板自由振動中的應用

伍洋 邢譽峰

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

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

中圖分類號： O323； O241.82 文獻標志碼：A文章編號1004-4523（2018）02-0343-09

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

引言

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

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

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

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

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

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

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

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

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

圖1Hermite插值基函數

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

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

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

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

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

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個頂點的直角坐標。此外，在薄板單元中需要計算到二階導數，因此需要用到如下鏈式法則：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形函數

與傳統等參單元類似，本文將在參考坐標系下構造單元的形函數，這樣形函數的參考域與幾何定義的參考域相同，最后在計算過程中通過參數變換可以得到物理域上的形函數。如圖2所示，單元在4個頂點各配置4個自由度，分別為w， wξ， wη， wξη， 而在各邊界內部每個節點配置兩個自由度w， wn， 即撓度和法向轉角， 各邊上的節點數目分別為M， N， P， Q。根據混合函數插值，那么可以得到每個自由度對應的形函數，如頂點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）其中上標V1代表頂點1，類似地可以得到其他頂點處的形函數，不再贅述。同理，也可以得到各邊界內部節點對應的形函數，如邊界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）

式中上標（1）代表邊界1，指標i代表邊（1）上的第i個節點。圖3給出了上述形函數的圖形。

由于薄板單元要滿足C1連續性，故單元內部升階譜形函數在單元邊界上必須滿足函數值以及導數值均為0。可以通過雅克比正交多項式的張量積形式來進行構造，即

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

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

式中J（4，4）m-1為雅克比正交多項式，權系數（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）形函數的正交性有利于提高單元結構矩陣的數值特性。其前兩個形函數如圖4所示， 可以看到，

圖3 邊界形函數

Fig.3 Edge shape functions圖4升階譜面函數

Fig.4Hierarchical face functions面函數在單元邊界上的函數值以及導數值均為0。

1.4數值離散

薄板的應變能與動能系數可以表示為U=D2S2wx22+2wy22+2ν2wx22wy2+

21-ν2wxy2dxdy

T0=12Sρhw2dxdy（15）式中ν為泊松比，D=Eh3/[12（1-ν2）]為薄板的彎曲剛度，E為彈性模量，ρ為密度，h為板的厚度。將上式離散即可得到單元剛度矩陣以及質量矩陣，為此，設單元的位移試函數為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ξη為頂點處的節點位移；w（j）i，w（j）n，j=1～4為邊內的節點位移；amn為對應升階譜面函數的廣義節點位移。定義如下節點位移向量= [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通過高斯-洛巴托積分，詳見文獻[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為沿著ξ，η方向的高斯-洛巴托積分系數；Gxx， Gyy， Gxy為微分矩陣。注意到是參考系下的節點位移向量，因而不能直接用于單元組裝與邊界條件的施加，為此需要將其轉化為x-y坐標系中的節點位移，如圖2（右）所示，為此定義如下節點位移向量= [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）同時可以得到可逆轉換關系=Q（21）將上式代入式（17）可以得到U=12TQ-THTBHQ-1

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

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

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

圖5扇形板

Fig.5A sectorial plate

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

圖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結論

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

參考文獻：

[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]. 北京：國防工業出版社， 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