曲梁结构由于其造型优美、施工方便、受力性能好等优点，被广泛运用于航天、机械、土木等领域。对曲梁结构的自由振动，国内外许多学者作了深入研究并提出若干分析方法，主要包括：有限元法[1—3]、有限差分法[4]、微分容积法[5]、动力刚度法[6—8]、伪谱法[9]等。这些方法中有限元法适用性最强，应用范围最广。
有限元求解精度取决于网格疏密和单元次数，网格越密，单元次数越高，精度越好，但所需计算量也随之快速增长，尤其是高阶频率和振型，为了达到满意的精度，往往需要花费很大的计算量[1]。因此，探索如何有效地提高有限元求解自由振动问题的精度和效率具有十分重要的意义，超收敛计算则是实现这一目标的有效途径之一。
目前结构自由振动问题的有限元超收敛计算大多从振型入手，先对振型进行超收敛修复，再将超收敛振型代入Rayleigh商得到频率的超收敛解。在这些超收敛计算方法中，比较典型的如在单元内进行二次划分并用内部结点位移进行插值校正的插值法[10]；选取若干单元构成单元片，在单元片上利用利用结点位移的超收敛性对振型进行最小二乘拟合，修复振型误差的SPRD法[11]；通过对应力进行插值，提高Rayleigh商中应变能计算精度从而提高频率精度的应力超收敛恢复法[12—13]；借助等几何分析方法构造高阶质量矩阵，提高频率计算精度的等几何分析法[14]；以及基于外推[15]和基于投影[16]的方法等。文[17]提出一种自由振动问题的p型超收敛算法，效果显著。
本文针对平面曲梁面内自由振动问题，推导了p型超收敛算法的算法公式，并进行了误差估计。本文工作是文[17]工作的延续，实现了p型超收敛算法由简单自由振动问题到平面曲梁这类更加复杂、更加通用问题的推广。数值算例显示，该法延续了在求解简单自由振动问题中的优秀表现，算法可靠、高效，能显著提高解答的精度、质量和收敛阶，再一次展示了p型超收敛算法在求解自由振动问题方面的优良特性。
1 曲梁自由振动分析模型
如图1所示平面曲梁，取轴线弧长坐标s，法向坐标n由轴线坐标s逆时针转90°得到，定义曲率中心在n的负方向时曲率为正。将曲梁截面形心处轴向、横向和转角位移分别记为u(s)、v(s)、θ(s)，截面内力N、Q、M以图1所示方向为正。
自由振动时，平面曲梁的几何方程、物理方程和平衡方程如下[9]：
图1 平面曲梁
Fig.1 Planar curved beam
式中：
EA为截面拉压刚度；EI为抗弯刚度；kGA为剪切刚度；R为曲率半径。
将几何方程和物理方程代入到平衡方程中，得到用位移表示的曲梁面内自由振动控制微分方程：
将式(4)写成向量形式：
式中：l为曲梁的轴线长度；L、m为相应的微分算子矩阵，其具体表达式可由式(4)得出，此处从略。由式(5)可以看出本质上曲梁面内自由振动是一个2阶向量型常微分方程广义特征值问题。该问题理论上有无穷多组解答
特征值按λn从小到大排序：
有限元求解时，先通过网格划分，对原结构进行离散，如图2所示，将结构划分为ne个单元。记
图2 网格划分
Fig.2 Finite element mesh
网格尺寸为：
式中：he为单元(e)的轴线长度，
对该问题，数值试验表明，有限元解的收敛阶取决于
中三个位移分量的形函数阶数的最低值，为使有限元自由度的求解效率得到充分发挥，应让三个位移分量的形函数取相同阶数。故本文有限元计算和超收敛计算中，三个位移分量形函数的阶数均取为相同。在任意单元(e)内，采用m次Hierarchical形函数将解答近似为：
式中：
为无量纲坐标；
分别为u、v、θ的内部自由度。
为n阶振型第e个结点的振型结点位移，将式(9)写成矩阵形式：
式中，
为单元(e)的自由度向量，具体为：
式中，
分别为u、v、θ的形函数矩阵，其具体表达式可由式(9)和式(11)得出，此处不再详细给出。
由变分原理[1]，有限元法将原问题离散为如下矩阵广义特征值问题：
式中：
为该网格下的整体刚度矩阵和质量矩阵，由单元刚度与质量矩阵集成获得；{Δh}为自由度向量，单元刚度与质量矩阵具体计算 如下：
求解该广义特征值问题，即可获得该网格下的有限元解
文[1]指出，对该向量型常微分方程特征值问题，频率和振型有如下收敛阶：
式中，m为单元次数。
文[18]中对该向量型问题对应的静力问题
只有一个分量的最简单情形，证明了结点位移具有h2m阶的超收敛性。数值算例表明，该结论对向量型特征值问题同样成立，即有：
振型结点位移的超收敛性是本文算法的基础。
2 超收敛求解思路
由式(5)可知，第n阶振型
在单元(e)上满足如下常微分方程边值问题，当网格足够密时，为该边值问题的唯一解：
由式(15)和式(17)知，振型结点位移和特征值均具有h2m阶超收敛性，在网格加密时，这些值将比内点位移更快地向精确解收敛，因此当网格足够密时，将有：
此时，单元(e)上的振型将近似于下列边值问题的解：
本文通过求解该线性边值问题获得单元(e)上振型的超收敛解，对该边值问题，本文直接采用高次单元进行求解。设单元次数为
单元上振型超收敛解为：
式中：
为单元上的
次形函数矩阵；下标b对应单元端部结点位移自由度；下标c对应单元内部自由度，各位移分量的具体表达式为：
该问题的有限元解方程为：
式中：
为单元动力刚度矩阵中端部自由度和内部自由度对应的子矩阵，单元动力刚度阵具体如下：
由式(25)解得：
从而得到单元(e)上振型的超收敛解：
对每个单元都进行上述超收敛修复，即可得到曲梁全域振型的超收敛解答。将振型超收敛解
代入到Rayleigh商，即可得到第n阶频率的超收敛解答：
3 误差估计
分别记振型有限元解和超收敛解的误差为：
其中：
单元上 d*是
的
次有限元解，由文献[1]的估计：
由式(18)和式(21)相减，得：
由式(15)、式(17)知：
故由式(35)知，当网格足够密时(密至每个单元的固端特征值均大于
时)，有：
从而有超收敛解的误差为：
而根据文献[1]，有限元解的误差：
易见，当m＞1的时候
d*收敛阶高于 dh，具有超收敛性。数值算例显示，将d*代入Rayleigh商得到的 λ*也有超收敛性。
4 数值算例
本文算法已编成Fortran程序，本节通过具体算例来展示本文方法的效果。为简单起见，计算中采用无量纲计算，为考察振型收敛阶，本文选取振型误差向量模为长度的最大模：
曲梁面内振动问题只有常截面圆弧曲梁才能推导出解析解，本文前两个算例选取不同边界条件下的常截面圆弧曲梁，用以展示算法的收敛阶，第三个算例是为了与文献中的计算结果进行对比，以展示算法的超收敛修复效果。
例1.常截面圆弧曲梁1
如图3所示半圆常截面圆弧曲梁，圆弧半径为R，曲梁为矩形截面，高h=R/3，宽b=h/2，截面形状系数k=5/6，弹性模量为E，剪切模量为
记无量纲频率
图3 例1常截面圆弧曲梁
Fig.3 Arch beam in Example 1
选取转角α为坐标变量(0≤α ≤π)，该圆弧曲梁的频率和振型可精确求出，参考文献[8]，其振型可假设为：
式中：n≥0为整数。将其代入式(4)中，可以得到关于
的齐次方程，由该齐次方程的系数矩阵奇异，可以推导出如下形式的频率方程：
借助卡丹公式，对每个n可以求得三个频率值。本结构前两阶频率为0，此处讨论第三阶频率，该阶频率对应n=2，频率方程见式(43)，频率解和振型解分别见式(44)和式(45)。
对该阶频率和振型分别用二次元(m=2)、三次元(m=3) 和四次元(m=4)进行计算，取
作振型归一化。考察曲梁自振频率误差
振型误差
和曲梁跨中截面结点C处(α=π/2)位移误差
的收敛阶，结果如图4和表1所示。结果表明频率和振型具有式(15)、式(16)所示的收敛阶，振型结点位移也具有式(17)所示的超收敛性。
图4 例1第三阶频率、第三阶振型和结点位移的收敛阶
Fig.4 Rate of convergence of the 3rd frequence, the 3rd mode and its nodal displacement in Example 1
为考察本文方法的精度和修复效果，将该曲梁均分为5个一次元(m=1)进行有限元求解，用二次元
进行超收敛修复，考察第三阶频率和振型，此时频率的有限元解为
(相对误差为154.8%)，超收敛解
(相对误差为12.9%)，图5给出了该阶振型的有限元解、超收敛解和精确解的比较，取
进行归一化。易见，本文方法对频率和振型的精度修复效果非常显著。
为考察本文方法的收敛阶，选取第三阶频率和振型，采用二次元(m=2)进行有限元计算，分别用三次元
、四次元
和五次元
进行超收敛修复，通过均匀二分加密网格统计频率和振型的收敛阶。计算结果如图6和表2所示。易见，相对于有限元解，超收敛解的精度和收敛阶有显著提升。本例超收敛计算中，由四次元
升到五次元
时，收敛阶保持不变，进一步研究发现，若继续提高
超收敛解的收敛阶仍与的收敛阶相同，这说明
的收敛阶与
收敛阶相同，进一步提高
超收敛解的收敛阶并不继续提高。
表1 例1第三阶振型及结点位移的收敛阶
Table 1 Rate of convergence of the 3rd frequency, the 3rd mode and its nodal displacement in Example 1
ne | |hΩ Ω- rcc || ||h-d d r || ||h-d d r 4 2.6205×10-1 — 2.8549×10-2 — 2.8549×10-2 8 2.0131×10-2 3.7 2.4594×10-3 3.5 3.7963×10-3 2.9 16 1.3422×10-3 3.9 1.6965×10-4 3.9 4.8478×10-4 3.0 32 8.5353×10-5 4.0 1.0885×10-5 4.0 6.0693×10-5 3.0 64 5.3582×10-6 4.0 6.8485×10-7 4.0 7.5874×10-6 3.0 80 2.1959×10-6 4.0 2.8074×10-7 4.0 3.8848×10-6 3.0 m=2 4 4 3 4 5.1139×10-3 — 9.1567×10-4 — 3.7427×10-3 8 9.2205×10-5 5.8 1.6207×10-5 5.8 2.1753×10-4 4.1 16 1.4947×10-6 5.9 2.6139×10-7 6.0 1.2742×10-5 4.1 32 2.3570×10-8 6.0 4.1168×10-9 6.0 7.7993×10-7 4.0 64 3.6837×10-10 6.0 6.4422×10-11 6.0 4.8475×10-8 4.0 80 9.5680×10-11 6.0 1.6815×10-11 6.0 1.9842×10-8 4.0 m=3 6 6 4 4 5.3399×10-5 — 7.4180×10-6 — 2.3241×10-4 8 2.2955×10-7 7.9 3.0684×10-8 7.9 7.4231×10-6 5.0 12 9.1140×10-9 8.0 1.2107×10-9 8.0 9.7187×10-7 5.0 16 9.2014×10-10 8.0 1.2170×10-10 8.0 2.2995×10-7 5.0 20 1.5393×10-10 8.0 2.0407×10-11 8.0 7.5232×10-8 5.0 24 3.6060×10-11 8.0 4.7521×10-12 8.0 3.0207×10-8 5.0 m=4 8 8 5
图5 例1第三阶振型有限元解与超收敛解的比较
Fig.5 The comparison of the FE solution with the recovered solution on the 3rd mode in Example 1
为将本文方法与SPRD法[11]进行对比，本文对第三阶频率和振型的二次元解用SPRD法进行超收敛修复，用相邻四个结点(单元两端结点外加离单元最近的两个结点)上的位移解作三次多项式插值得单元上振型的超收敛解，将振型超收敛解代入Rayleigh商得频率超收敛解。具体结果如表3所示。易见，SPRD法三次元修复解与本文三次元
修复解具有相同的收敛阶，只是本文方法的解答精度比SPRD法稍高一些。
图6 例1第三阶特征对有限元解与超收敛解收敛阶
Fig.6 Rate of convergence of the FE solution and recovered solutions on the 3rd eigenpair in Example 1
表2 例1第三阶特征对有限元解与超收敛解收敛阶
Table 2 Rate of convergence of the FE solution and recovered solutions on the 3rd eigenpair in Example 1
2 m= ne | |h Ω Ω- r || ||h-d d r 4 2.6205×10-1 — 2.8549×10-2 — 8 2.0131×10-2 3.7 3.7963×10-3 2.9 16 1.3422×10-3 3.9 4.8478×10-4 3.0 32 8.5353×10-5 4.0 6.0693×10-5 3.0 64 5.3582×10-6 4.0 7.5874×10-6 3.0 80 2.1959×10-6 4.0 3.8848×10-6 3.0理论值 4 3 3 m= ne *| |Ω Ω- r *|| ||-d d r 4 1.6240×10-2 — 2.8549×10-2 — 8 2.9486×10-4 5.8 2.4594×10-3 3.5 16 2.8338×10-6 6.7 1.6965×10-4 3.9 32 2.9635×10-8 6.6 1.0970×10-5 4.0 64 3.9263×10-10 6.2 6.9733×10-7 4.0 80 1.0010×10-10 6.1 2.8630×10-7 4.0 理论值 6 4 4 m= ne *| |Ω Ω- r *|| ||-d d r 4 1.1863×10-2 — 2.8549×10-2 — 8 2.0486×10-4 5.9 2.4594×10-3 3.5 16 1.3430×10-6 7.3 1.6965×10-4 3.9 32 6.0706×10-9 7.8 1.0885×10-5 4.0 64 2.3540×10-11 8.0 6.8485×10-7 4.0 80 3.9800×10-12 8.0 2.8074×10-7 4.0理论值 8 4 5 m= ne *| |Ω Ω- r *|| ||-d d r 4 1.1815×10-2 — 2.8549×10-2 — 8 2.0463×10-4 5.9 2.4594×10-3 3.5 16 1.3420×10-6 7.3 1.6965×10-4 3.9 32 6.0670×10-9 7.8 1.0885×10-5 4.0 64 2.3530×10-11 8.0 6.8485×10-7 4.0 80 3.9800×10-12 8.0 2.8074×10-7 4.0理论值 8 4
例2.常截面圆弧曲梁2
如图7所示，本例亦为半圆常截面圆弧曲梁，所有参数均与例1相同，只是两端的支座变成了水平滑动支座。
该曲梁振型可写为下列形式：
表3 例1第三阶特征对SPRD法收敛阶
Table 3 Rate of convergence of the recovered solutions on the 3rd eigenpair in Example 1 with SPRD method
?
图7 例2常截面圆弧曲梁
Fig.7 Arch beam in Example 2
式中，n≥0为整数。
该曲梁第一阶频率为Ω1=0，此时的振型即水平方向的刚体平动(具体为n=1，u=sin(α)、
第四阶频率
振型为n=0，u=0、v=1、θ=0，即整个半圆弧发生沿着径向的伸缩，有限元多项式形函数能准确描述该振型，因此该阶振型的有限元解为精确解。
将曲梁均分为32个单元(ne=32)，采用一次元(m=1)进行有限元计算，用二次元
进行修复。表4列出该网格下前12阶频率和振型的有限元解与超收敛解的精度对比。易见，对一次元， 虽然振型结点位移不具有超收敛性，本文方法仍能显著提高解答(尤其是频率)的精度。
本例选取该曲梁的高阶频率来研究算法收敛阶，选取第27阶频率和振型，此时n=15，用例1中类似方法可计算精确解，该阶频率和振型的精确解为：
取三次元(m=3)进行有限元计算，分别用四次元
、五次元
和六次元
来进行超收敛修复，取v(α=0)=1进行归一化，计算结果如图8和表5所示。结果表明本文算法对高阶频率和振型仍能够有效修复。
例3.一般轴线曲梁
文[3]和文[19]分别研究了图9所示的抛物线、三角正弦线和椭圆线三类轴线曲梁的前四阶频率，其轴线方程分别如下：
其中，
定义高跨比为长细比
为
截面剪切刚度与轴向刚度之比为
本例无量纲频率记为
表4 例2前12阶特征对有限元解与超收敛解误差
Table 4 Error of the FE solution and recovered solutions on the first 12 eigenpairs in Example 2
阶数 Ω h Ω | |h Ω Ω- || ||h-d d *Ω *|- |ΩΩ*| |(%)| |hΩ ΩΩ Ω--*|| ||-d d *|| ||(%)|| ||h--d d d d n 1 0.0000000000 0.2406633400 2.41×10-1 1.47×10-3 0.0036229587 3.62×10-3 1.5 2.79×10-4 18.9 1 2 2.5023890834 2.6333380484 1.31×10-1 4.81×10-3 2.5025246366 1.36×10-4 0.1 8.71×10-4 18.1 2 3 6.5859454955 6.8160091249 2.30×10-1 1.08×10-2 6.5864362731 4.91×10-4 0.2 1.53×10-3 14.1 3 4 10.392304845 10.392304845 4.99×10-13 3.48×10-14 10.392304845 0.00×100 — 3.48×10-14 — 0 5 11.633033136 12.035028151 4.02×10-1 2.00×10-2 11.634507076 1.47×10-3 0.4 2.83×10-3 14.1 4 6 14.628136082 14.630131136 2.00×10-3 1.52×10-3 14.628136272 1.91×10-7 0.0 2.70×10-4 17.8 1 7 17.270129270 17.928031623 6.58×10-1 3.54×10-2 17.273953249 3.82×10-3 0.6 5.63×10-3 15.9 5 8 23.169351444 23.195502842 2.62×10-2 6.17×10-3 23.169363367 1.19×10-5 0.0 8.40×10-4 13.6 2 9 23.263957307 24.275198368 1.01×100 4.44×10-2 23.272759951 8.80×10-3 0.9 9.78×10-3 22.0 6 10 29.470096838 30.945332909 1.48×100 5.99×10-2 29.488468798 1.84×10-2 1.2 1.64×10-2 27.3 7 11 32.809197811 32.909640208 1.00×10-1 1.40×10-2 32.809322307 1.24×10-4 0.1 1.54×10-3 11.0 3 12 35.799123270 37.861888913 2.06×100 7.44×10-2 35.834527464 3.54×10-2 1.7 2.37×10-2 31.9 8
图8 例2第27阶特征对有限元解与超收敛解收敛阶
Fig.8 Rate of convergence of the FE solution and recovered solutions on the 27th eigenpair in Example 2
表5 例2第27阶特征对有限元解与超收敛解收敛阶
Table 5 Rate of convergence of the FE solution and recovered solutions on the 27th eigenpair in Example 2
3 m= ne | |h Ω Ω- r || ||h-d d r 8 3.5762×100 9.9343×100 16 2.7289×10-1 3.7 2.2938×10-2 8.8 32 5.3927×10-3 5.7 2.5326×10-3 3.2 64 8.9270×10-5 5.9 1.5439×10-4 4.0 80 2.3563×10-5 6.0 6.2989×10-5 4.0 100 6.2043×10-6 6.0 2.5759×10-5 4.0 理论值 6 4 4 m= ne *| |Ω Ω- r *|| ||-d d r 8 3.0414×100 1.0961×101 16 1.3966×10-2 7.8 2.2938×10-2 8.9 32 4.7421×10-5 8.2 2.0408×10-4 6.8 64 1.9294×10-7 7.9 5.6747×10-6 5.2 80 3.2539×10-8 8.0 1.8226×10-6 5.1 100 5.4775×10-9 8.0 5.9733×10-7 5.0 理论值 8 5
续表
5 m= ne *| |Ω Ω- r *|| ||-d d r 8 3.0789×100 1.1835×101 16 4.1780×10-3 9.5 2.2938×10-2 9.0 32 4.2711×10-7 13.3 7.2806×10-5 8.3 64 2.8419×10-10 10.6 1.0749×10-6 6.1 80 2.9800×10-11 10.1 2.8069×10-7 6.0 100 3.1974×10-12 10.0 7.4345×10-8 6.0 理论值 10 6 6 m= ne *| |Ω Ω- r *|| ||-d d r 8 3.0795×100 1.1889×101 16 3.9569×10-3 9.6 2.2938×10-2 9.0 32 1.6648×10-7 14.5 7.1572×10-5 8.3 64 1.9597×10-11 13.1 1.0003×10-6 6.2 80 1.1612×10-12 12.7 2.5941×10-7 6.0 100 7.2109×10-14 12.5 6.7551×10-8 6.0 理论值 12 6
图9 三类轴线曲梁
Fig.9 Three types of curved beams
本文对这三类曲梁计算其前4阶频率，沿x轴均匀划分为32个单元，分别采用一次元
进行有限元计算，后采用二次元
进行超收敛修复和采用二次元
进行有限元计算，后采用三次元
进行超收敛修复，计算结果如表6所示。易见，在该网格下，若采用一次元进行有限元计算，部分频率仍有较大误差，但采用二次元进行修复后与文献[3]的结果相比误差均在0.2%以下，其精度与直接采用二次元
进行有限元计算的结果精度较为接近。而在二次元有限元解答的基础上使用三次元进行修复，得到的频率解答与文 献[3]几乎完全一致。表7给出了上述32个单元网格下三种轴线曲梁前40阶振型有限元计算和超收敛计算的总耗时。易见，由于本文的超收敛计算是在后处理阶段的一系列单元局部小规模问题的线性计算，耗时非常少，但解答精度得到大幅提升，一次元求解基础上用二次元作超收敛恢复，得到的解答结果与直接用二次元求解的精度较为接近，但计算时间不到二次元求解的二分之一，充分体现出本文方法对精度恢复的高效特性。
表6 例3三类轴线曲梁的无量纲频率
Table 6 Non-dimensional frequencies for three curved beams of Example 3
曲梁 阶数 文献[3]文献[19]32 1 2 ne m m ===，， 32 2 3 ne m m ===，， 有限元解 误差/(%) 超收敛解 误差/(%) 有限元解 误差/(%) 超收敛解 抛物线 铰支-铰支 f=0.3 SR=75 kG/E=0.31 21.759 21.83 25.3654 16.57 21.7755 0.08 21.7612 0.01 21.7590 2 55.493 56 62.7888 13.15 55.5478 0.10 55.5040 0.02 55.4931 3 100.701 102.3 113.5027 12.71 100.8287 0.13 100.7365 0.04 100.7017 4 113.302 113.4 113.6644 0.32 113.3626 0.05 113.3032 0.00 113.3018 三角函数线 固支-固支 ε=0.5, f=0.1 SR=100 kG/E=0.31 56.083 56.3 63.4051 13.06 56.0998 0.03 56.0964 0.02 56.0829 2 66.047 66.14 67.4952 2.19 66.1787 0.20 66.0496 0.00 66.0473 3 113.406 114.3 127.4757 12.41 113.3699 -0.03 113.4482 0.04 113.4063 4 179.264 181.7 203.8437 13.71 179.446 0.10 179.3757 0.06 179.2652 椭圆线 铰支-固支 ε=0.5, f=0.2 SR=50 kG/E=0.31 34.892 35.25 36.6147 4.94 34.9018 0.03 34.8948 0.01 34.8924 2 56.766 57.11 57.4903 1.28 56.8260 0.11 56.7670 0.00 56.7655 3 81.42 83 84.6606 3.98 81.3835 -0.04 81.4277 0.01 81.4203 4 124.288 128.2 130.0805 4.66 124.5121 0.18 124.3096 0.02 124.2877
表7 例3前40阶振型计算总耗时 /s
Table 7 Total computation time for the first 40 modes of
Example 3
有限元解1 (m=1)超收敛解1 (2 m= )有限元解2 (m=2)超收敛解2 (3 m= )a 6.112 0.028 13.135 0.056 b 6.281 0.021 13.266 0.053 c 6.249 0.019 13.044 0.053
以上算例充分表明本文方法稳定、高效，能够在有限元解的基础上，显著提高解答的精度、质量和收敛阶。数值算例表明，本文超收敛解答具有如表8所示的收敛规律。
表8 超收敛解的收敛阶
Table 8 The convergence order of the recovered solutions
*ω *d r min(2,4)m m min( 1,2 )m m+
5 结论
本文的超收敛算法基于有限元解答中频率和振型结点位移固有的的超收敛性质，能够在有限元解的基础上有效地提高解答的精度和质量，在较粗的网格下即可获得较好精度的频率和振型，思路巧妙明晰，颇具潜力，值得进一步推广。
[1]Strang G, Fix G J.An analysis of the finite element method [M].Prentice-Hall: Englewood Cliffs, NJ, 1973.
[2]Tüfekçi E, Arpaci.Exact solution of in-plane vibrations of circular arches with account taken of axial extension, transverse shear and rotatory inertia effects [J].Journal of Sounds & Vibration, 1998, 209(5): 845—856.
[3]Yang F, Sedaghati R, Esmailzadeh E.Free in-plane vibration of general curved beams using finite element method [J].Journal of Sound & Vibration, 2008, 318(4): 850—867.
[4]Ball R.Dynamic analysis of rings by finite differences [J].Journal of the Engineering Mechanics Division, 1967, 93(1): 1—10.
[5]Kang K, Bert C W, Striz A G.Vibration analysis of shear deformable circular arches by the differential quadrature method [J].Journal of Sound & Vibration, 1995, 183(2): 353—360.
[6]Tseng T P, Huang C S, Lin C J.Dynamic stiffness analysis for in-plane vibration of arches with variable curvature [J].Journal of Sound & Vibration, 1997, 207(1): 15—31.
[7]叶康生, 赵雪健.动力刚度法求解平面曲梁面外自由振动问题[J].工程力学, 2012, 29(3): 1—8.Ye Kangsheng, Zhao Xuejian.Dynamic stiffness method for out-of-plane free vibration analysis of planar curved beams [J].Engineering Mechanics, 2012, 29(3): 1—8.(in Chinese)
[8]Howson W P, Jemah A K.Exact dynamic stiffness method for planar natural frequencies of curved Timoshenko beams [J].Journal of Mechanical Engineering Science, 1999, 213(7): 687—696.
[9]Lee J.In-plane free vibration analysis of curved Timoshenko beams by the pseudospectral method [J].Ksme International Journal, 2003, 17(8): 1156—1163.
[10]林群, 杨一都.解二阶椭圆本征值问题的有限元插值校正方法[J].计算数学, 1992, 14(3): 334—338.Linqun, Yang Yidu.The finite element interpolated correction method for second order elliptic eigenvalue problems [J].Mathematica Numerica Sinica, 1992, 14(3): 334—338.(in Chinese)
[11]Wiberg N E, Bausys R, Hager P.Improved eigenfrequencies and eigenmodes in free vibration analysis [J].Computers & Structures, 1999, 73(1/2/ 3/4/5): 79—89.
[12]龚国庆, 范子杰, 刘寒冰.基于应力超收敛恢复技术的广义特征值问题后验误差估计[J].工程力学, 2004, 21(3): 111—117.Gong Guoqing, Fan Zijie, Liu Hanbing.A posteriori error estimation based on stress super-convergence recovery technique for generalized eigenvalue problems [J].Engineering Mechanics, 2004, 21(3): 111—117.(in Chinese)
[13]Avrashi J, Cook R D.New error estimation for C0 eigenproblems in finite element analysis [J].Engineering Computations, 1993, 10(3): 243—256.
[14]Wang D, Liu W, Zhang H.Novel higher order mass matrices for isogeometric structural vibration analysis [J].Computer Methods in Applied Mechanics & Engineering, 2013, 260(9): 92—108.
[15]Li M X, Lin Q, Zhang S H.Extrapolation and superconvergence of the Steklov eigenvalue problem [J].Advances in computational mathematics, 2010, 33(1): 25—44.
[16]Huang P Z.Superconvergence of a stabilized approximation for the stokes eigenvalue problem by projection method [J].Applications of mathematics, 2014, 59(4): 361—370.
[17]叶康生, 曾强.结构自由振动问题有限元新型超收敛算法研究[J].工程力学, 2017, 34(1): 45—50, 68.Ye Kangsheng Zeng qiang.A new superconvergent recovery method for FE analysis on structural free vibration problems [J].Engineering Mechanics, 2017, 34(1): 45—50, 68.(in Chinese)
[18]Douglas J, Dupont T.Galerkin approximations for the two point boundary problems using continuous piecewise polynomial spaces [J].Numerical Mathematics, 1974, 22(2): 99—109.
[19]Oh S J, Lee B K, Lee I W.Natural frequencies of non-circular arches with rotatory inertia and shear deformation [J].Journal of Sound & Vibration, 1999, 219(1): 23—33.