平面曲梁结构在现代土木、航天、机械等领域具有广泛的应用，其自由振动是结构分析的重要内容。该问题目前的分析方法主要有：有限元法^[1-3]、有限差分法^[4]、微分容积法^[5-6]、伪谱法^[7]、动力刚度法^[8-10]等，其中有限元法应用最为广泛。

有限元法求解自由振动问题时，计算精度由单元次数和网格划分决定。网格越密，单元次数越高，计算精度越好，但计算量亦随之大幅攀升，高阶频率和振型尤甚^[1]。因此，研究如何有效提高有限元求解自由振动问题的精度和效率具有十分重要的意义，超收敛计算则是解决这一难题的有效途径。

目前结构自由振动问题的各类有限元超收敛求解方法，基本都是先对振型进行超收敛修复，再将修复后的振型代入Rayleigh商得到频率的超收敛解，区别主要在于对振型的超收敛修复做法不同。这些方法中最典型的是文献[11]提出的SPRD法，该法利用振型结点位移的超收敛性，在多个单元组成的小片上通过对振型用更高次多项式对结点位移进行多项式拟合获得振型的超收敛解。文献[12]的做法与文献[11]做法类似，亦是基于振型结点位移的超收敛性，在几个单元构成的小片上通过对振型用更高次多项式对结点位移作多项式插值获得振型的超收敛解；文献[13-14]仅修复振型所对应的应力，并代入Rayleigh商的应变能部分，通过应变能计算精度的提高获得频率的超收敛解；文献[15-16]利用等几何分析法，通过构造高阶质量矩阵，获得频率的高精度解；另外还有文献[17]基于外推的方法和文献[18]基于投影的方法。

文献[19]对自由振动问题提出一种p型超收敛算法，结果显示，该法简单、高效，能显著提升频率和振型的精度和收敛阶。文献[20]将该法应用于平面曲梁的静力分析中，文献[21]将该法延拓至平面曲梁的面内自由振动分析中，文献[22]将该法延拓至欧拉梁的弹性稳定分析中，文献[23]将该法延拓至非线性常微分方程边值问题的分析中。本文的工作是文献[19-23]工作的延续，是文献[21]工作的后半部分，进一步将p型超收敛算法应用于平面曲梁的面外自由振动分析中。数值算例表明，该法仍延续了求解前述问题时的优秀表现，算法稳定、高效，能显著提高解答的精度、质量和收敛阶，再次展现了p型超收敛算法在求解自由振动问题时的优良特性。

1   平面曲梁自由振动分析模型

平面曲梁的轴线为平面曲线，所有截面关于轴线所在平面对称。如图1所示，建立曲梁分析的$xy{\textit{z}}$坐标系，$x$、$y$分别沿轴线切向、法向，${\textit{z}}$垂直于轴线所在平面。记$s$为轴线弧长坐标。曲梁面外自由振动时，分别记面外位移、$y$轴转角位移和扭转角位移为$w (s)$、$\theta (s)$、$\phi (s)$，对应内力$Q$、$M$、$T$如图1所示。

图  1  平面曲梁

Figure  1.  Planar curved beam

下载: 全尺寸图片 幻灯片

面外自由振动时，曲梁的几何方程、物理方程和平衡方程分别为^[24]：

式中：$(\;)' = {{{\rm{d}}(\;)} / {{\rm{d}}s}}$；$E$为弹性模量；$G$为剪切模量；$\rho$为材料密度；$A$为截面面积；${I_y}$为截面关于$y$轴的惯性矩；${I_{\rm p}}$为截面极惯性矩；k为截面形状系数；$J$为截面扭转常数；R为轴线曲率半径，当轴线曲率中心位于$y$轴正向时R为正，位于$y$轴负向时R为负。

将物理方程和几何方程代入平衡方程，得到用位移表示的平面曲梁面外自由振动的控制微分方程为：

将其写为向量形式：

式中：$l$为曲梁的轴线长度；L、m为相应的微分算子矩阵，其具体表达式可由式(4)导出，此处不再详细给出。该问题理论上有无穷多组解答$({\lambda _n},\;{{{d}}_n}(s)),$$\;n = 1,2, \cdots$，特征值${\lambda _n}$按从小到大排序：

用有限元法求解该问题时，先对结构进行有限元离散，将结构划分为$ne$个单元，如图2所示。

图  2  有限元网格划分

Figure  2.  Finite element mesh

下载: 全尺寸图片 幻灯片

记网格尺寸：

式中，${h_e}$为单元$(e)$的轴线长度，${h_e} = {s_{e + 1}} - {s_e}$。虽然理论上${{d}} = {\{ w,\theta ,\phi \} ^{\rm{T}}}$中三个位移分量的形函数阶数可取为不同，但数值结果表明，有限元解的收敛阶取决于这三个位移分量形函数阶数中的最低值，为使有限元自由度的求解效率得到充分发挥，应让三个位移分量的形函数阶数取为相同，故本文有限元计算和超收敛计算中三个位移分量的形函数阶数均取为相同。在任意单元$(e)$内，采用$m$次Hierarchical形函数将解答近似为：

式中：$\xi = ({{s - {s_e})} / {{h_e}}}$，为单元$(e)$上的无量纲坐标；${w_e}$、${\theta _e}$、${\phi _e}$为第e个结点上的结点位移；${a_i}$、${b_i}$、${c_i}$$(i = 1,2, \cdots ,m - 1)$分别为$w$、$\theta$、$\phi$的单元内部自由度。将式(9)写成矩阵形式：

式中，${ { \varDelta} ^e}$为单元$(e)$的自由度向量，具体为：

${{{ N}_w}}$、${{{ N}_\theta }}$、${{{ N}_\phi }}$分别为$w$、$\theta$、$\phi$的形函数矩阵，其具体表达式可由式(9)和式(11)导出。

由变分原理^[1]，有限元法将原问题离散为如下矩阵广义特征值问题：

式中：${{ K}^h}$和${{ M}^h}$分别为该网格下的整体刚度矩阵和整体质量矩阵，分别由单元刚度矩阵和单元质量矩阵集成获得；$\ {{\varDelta }}^h\ $为该网格下的结构整体自由度向量。单元刚度矩阵和单元质量矩阵具体计算如下：

求解上述矩阵广义特征值问题即可获得该网格下的有限元解$(\lambda _n^h,\;{{d}}_n^h(s))$，$n = 1,2, \cdots$。

对该向量型常微分方程特征值问题，文献[1]估计有限元解中频率和振型分别具有如下收敛阶：

式中，$m$为单元多项式次数。

文献[25]对该向量型问题对应的静力问题${{Ld}} = {{f}}$只有一个分量的最简单情形，证明当采用$m$次元求解时，结点位移具有${h^{2m}}$阶的超收敛性。数值算例表明，该结论对特征值问题同样成立，在自由振动问题中，振型结点位移同样具有${h^{2m}}$阶的超收敛性，即：

振型结点位移的超收敛性是本文算法的基础。

2   超收敛求解思路

由式(5)可知，在单元$(e)$上，第$n$阶振型${{{d}}_n}(s)$满足如下常微分方程边值问题：

且当网格足够密时，为该边值问题的唯一解。

由式(15)和式(17)知，频率和振型结点位移均具有${h^{2m}}$阶的超收敛性。故在网格加密过程中，这些值将比内点位移更快地向精确解收敛，当网格足够密时，有：

此时，单元$(e)$上的振型将近似于下列线性边值问题的解：

本文通过求解该线性边值问题获得单元$(e)$上振型的超收敛解。对该线性边值问题，本文直接采用一个高次单元进行求解。设高次单元次数为$\bar m(\bar m > m)$，单元$(e)$上振型的超收敛解为：

式中：${\bar { N}^e}$为单元上$\bar m$次形函数矩阵；下标b对应单元端部结点位移自由度；下标c对应单元内部自由度。三个位移分量的表达式为：

该线性边值问题的有限元方程为：

式中：${{\bar { k}}_{cb}^{*e}}$、${{\bar { k}}_{cc}^{*e}}$为单元动力刚度矩阵中端部自由度和内部自由度对应的子矩阵。单元动力刚度阵为：

由式(25)解得：

将式(27)代入式(22)得单元$(e)$上振型的超收敛解为：

采用上述超收敛过程修复所有单元，即可获得第$n$阶振型的超收敛解。将振型的超收敛解${{{d}}^*} = {{\rm{\{ }}{w^*},{\theta ^*},{\phi ^*}{\rm{\} }^{\rm {T}}}}$代入Rayleigh商，则可得频率的超收敛解为：

3   误差估计

设振型有限元解的误差为：

振型超收敛解的误差为：

式中：

由于${{{d}}^*}$是${{\bar d}}$在单个单元上的$\bar m$次有限元解，由文献[1]的估计：

将式(18)和式(21)相减得：

由式(15)和式(17)知：

由式(35)知，若网格划分足够密时(密至每个单元的固端特征值均大于$\lambda _n^h$)，有：

从而有：

而振型有限元解的收敛阶为：

易见，当$m > 1$时，$\min (\bar m + 1,2m) > m + 1$，振型超收敛解${{{d}}^*}$的收敛阶高于有限元解${{{d}}^h}$，具有超收敛性。数值算例表明，将${{{d}}^*}$代入Rayleigh商后得到的频率解${\lambda ^*}$亦具有超收敛性。

4   数值算例

本文算法已编成Fortran程序。下面通过3个数值算例来展示本文算法的计算精度和修复效果。为简单起见，三个算例均采用无量纲计算。为考察振型收敛阶，本文选取振型误差向量模为长度的最大模：

平面曲梁面外自由振动问题只有常截面圆弧曲梁具有解析解。为展现本文算法对收敛阶的提升，例1和例2选取了不同边界条件的两个常截面圆弧曲梁；为展现本文算法的修复效果，例3与已有文献的计算结果进行了对比。

例1. 常截面圆弧曲梁1

如图3(a)所示轴线为半个圆周的常截面圆弧曲梁，圆弧半径为$R$，截面形状为矩形，如图3(b)所示。高$t = R/5$，宽$b = t/2$，截面形状系数$k = 5/6$，弹性模量为$E$，泊松比$\nu = 0.3$。为计算简便，采用无量纲频率$\varOmega = \omega {R^2}\sqrt {\rho A/EI}$。

图  3  例1常截面圆弧曲梁

Figure  3.  Uniform arch beam in Example 1

下载: 全尺寸图片 幻灯片

坐标参数选为转角$\alpha$$(0 \leqslant \alpha \leqslant \pi )$，边界条件为：

可求出该圆弧曲梁频率和振型的精确解，根据文献[26]，设其振型为：

式中，$n \geqslant 0$且为整数。将式(42)代入式(4)，化简可得关于${\{ {w_0},{\theta _0},{\phi _0}\} ^{\rm T}}$的齐次方程组。由于振型存在非零解，故该齐次方程组的系数矩阵奇异，其行列式值为零，从而可推出如下的频率方程：

由该频率方程可求出每个$n$值对应的三个频率值。本例第一阶、第二阶频率均为0，故对第三阶频率求解。第三阶频率对应$n = 2$，频率方程如式(44)所示，求解该方程可得如式(45)所示频率。将求出的频率代入${\{ {w_0},{\theta _0},{\phi _0}\} ^{\rm T}}$的齐次方程组，可求得对应的振型，如式(46)所示，该阶振型取$w(\alpha = 0{\rm{)}} = {\rm{1}}$进行归一化。

首先考察频率、振型和振型结点位移有限元解的超收敛性。对该阶振型分别采用二次元$(m = 2)$、三次元$(m = 3)$和四次元$(m = 4)$进行有限元求解，计算频率误差$|\varOmega {\rm{ - }}{\varOmega ^h}|$，振型误差$||{{d}}{\rm{ - }}{{{d}}^h}||$和跨中截面结点$C$处($\alpha = \pi /2$)位移误差$||{{{d}}_c}{\rm{ - }}{{d}}_c^h||$并进行收敛阶的统计，结果如图4和表1所示。结果表明频率和振型结点位移具有$2m$阶的超收敛性，而振型的收敛阶仅为$m + 1$阶。

图  4  例1第三阶频率、振型和结点位移的收敛阶

Figure  4.  Rate of convergence of 3rd frequence, mode and its nodal displacement in Example 1

下载: 全尺寸图片 幻灯片

表  1  例1第三阶频率、振型和结点位移的收敛阶

Table  1.  Rate of convergence of 3rd frequency, mode and its nodal displacement in Example 1

ne	$|\varOmega {\rm{ - }}{\varOmega ^h}|$	r	$||{{{d}}_c}{\rm{ - }}{{d}}_c^h||$	r	$||{{d}}{\rm{ - }}{{{d}}^h}||$	r
16	2.8821×10−3	−	7.5917×10−6	−	4.9237×10−4	−
32	1.8948×10−4	3.9	5.7596×10−7	3.7	6.2308×10−5	3.0
64	1.1998×10−5	4.0	3.7692×10−8	3.9	7.8069×10−6	3.0
128	7.5231×10−7	4.0	2.3827×10−9	4.0	9.7640×10−7	3.0
256	4.7058×10−8	4.0	1.4934×10−10	4.0	1.2207×10−7	3.0
512	2.9417×10−9	4.0	9.3404×10−12	4.0	1.5259×10−8	3.0
m=2		4		4		3
ne	$|\varOmega {\rm{ - }}{\varOmega ^h}|$	r	$||{{{d}}_c}{\rm{ - }}{{d}}_c^h||$	r	$||{{d}}{\rm{ - }}{{{d}}^h}||$	r
16	3.2911×10−6	−	2.0671×10−7	−	1.3848×10−5	−
32	5.2669×10−8	6.0	3.3117×10−9	6.0	8.1494×10−7	4.1
64	8.2793×10−10	6.0	5.2071×10−11	6.0	5.0101×10−8	4.0
128	1.2956×10−11	6.0	8.1489×10−13	6.0	3.1181×10−9	4.0
256	2.0251×10−13	6.0	1.2738×10−14	6.0	1.9468×10−10	4.0
512	3.1645×10−15	6.0	1.9905×10−16	6.0	1.2164×10−11	4.0
m=3		6		6		4
ne	$|\varOmega {\rm{ - }}{\varOmega ^h}|$	r	$||{{{d}}_c}{\rm{ - }}{{d}}_c^h||$	r	$||{{d}}{\rm{ - }}{{{d}}^h}||$	r
16	2.0402×10−9	−	1.2643×10−10	−	2.3812×10−7	−
32	8.0841×10−12	8.0	4.9957×10−13	8.0	7.3852×10−9	5.0
64	3.1692×10−14	8.0	1.9571×10−15	8.0	2.3031×10−10	5.0
128	1.2390×10−16	8.0	7.6504×10−18	8.0	7.1933×10−12	5.0
256	4.8646×10−19	8.0	2.9890×10−20	8.0	2.2476×10−13	5.0
512	1.8912×10−21	8.0	1.1676×10−22	8.0	7.0236×10−15	5.0
m=4		8		8		5

下载: 导出CSV 

| 显示表格

为考察p型超收敛算法的计算精度和修复效果，将曲梁均匀划分为5个单元，先用一次元$(m = 1)$进行有限元求解，再用二次元$(\bar m = 2)$进行超收敛修复。得到第三阶频率的有限元解为$\varOmega _3^h =$${\rm{8}}{\rm{.398}}$(相对误差为300.5%)，超收敛解为$\varOmega _3^{\rm{*}} = {\rm{2}}{\rm{.664}}$(相对误差为27.0%)。图5为第三阶振型的有限元解、超收敛解和精确解的对比。可见，本文算法可以显著提升频率和振型的计算精度。

图  5  例1第三阶振型的有限元解与超收敛解比较

Figure  5.  Comparison of FE solution with recovered solution on 3rd mode in Example 1

下载: 全尺寸图片 幻灯片

为考察p型超收敛算法的收敛阶，对第三阶频率和振型，先用二次元$(m = 2)$进行有限元计算，再用三次元$(\bar m = 3)$、四次元$(\bar m = 4)$和五次元$(\bar m = 5)$进行超收敛修复并统计收敛阶，计算结果如图6和表2所示。易见，超收敛修复后，精度和收敛阶与有限元解相比均有显著提升。同时，由表2可以看出，当从四次元$(\bar m = 4)$升至五次元$(\bar m = 5)$时，频率和振型的收敛阶未随之提升。理论分析和数值算例均表明，当$\bar m > 2m$时，超收敛解的收敛阶与$\bar m = 2m$时相同，不会随$\bar m$的进一步提高而增加。

图  6  例1第三阶特征对有限元解与超收敛解收敛阶

Figure  6.  Rate of convergence of FE solution and recovered solutions on 3rd eigenpair in Example 1

下载: 全尺寸图片 幻灯片

表  2  例1第三阶特征对有限元解与超收敛解收敛阶

Table  2.  Rate of convergence of FE solution and recovered solutions on 3rd eigenpair in Example 1

$m = 2$
$ne$	$|\varOmega {\rm{ - }}{\varOmega ^h}|$	$r$	$||{{d}} - {{{d}}^h}||$	$r$
16	2.8821×10−3	−	4.9237×10−4	−
32	1.8948×10−4	3.9	6.2308×10−5	3.0
64	1.1998×10−5	4.0	7.8069×10−6	3.0
128	7.5231×10−7	4.0	9.7640×10−7	3.0
256	4.7058×10−8	4.0	1.2207×10−7	3.0
512	2.9417×10−9	4.0	1.5259×10−8	3.0
理论值		4		3
$\bar m = 3$
$ne$	$|\varOmega {\rm{ - }}{\varOmega ^{\rm{*}}}|$	$r$	$||{{d}} - {{{d}}^{\rm{*}}}||$	$r$
16	4.1065×10−5	−	1.4864×10−4	−
32	2.5175×10−7	7.3	9.6899×10−6	3.9
64	1.6723×10−9	7.2	6.1437×10−7	4.0
128	1.6325×10−11	6.7	3.8746×10−8	4.0
256	2.1574×10−13	6.2	2.4284×10−9	4.0
512	3.2163×10−15	6.1	1.5188×10−10	4.0
理论值		6		4
$\bar m = 4$
$ne$	$|\varOmega {\rm{ - }}{\varOmega ^{\rm{*}}}|$	$r$	$||{{d}} - {{{d}}^{\rm{*}}}||$	$r$
16	3.7774×10−5	−	1.4864×10−4	−
32	1.9909×10−7	7.6	9.6899×10−6	3.9
64	8.4441×10−10	7.9	6.1233×10−7	4.0
128	3.3689×10−12	8.0	3.8378×10−8	4.0
256	1.3230×10−14	8.0	2.4003×10−9	4.0
512	5.1747×10−17	8.0	1.5004×10−10	4.0
理论值		8		4
$\bar m = 5$
$ne$	$|\varOmega {\rm{ - }}{\varOmega ^{\rm{*}}}|$	$r$	$||{{d}} - {{{d}}^{\rm{*}}}||$	$r$
16	3.7772×10−5	−	1.4864×10−4	−
32	1.9908×10−7	7.6	9.6899×10−6	3.9
64	8.4438×10−10	7.9	6.1233×10−7	4.0
128	3.3687×10−12	8.0	3.8378×10−8	4.0
256	1.3229×10−14	8.0	2.4003×10−9	4.0
512	5.1745×10−17	8.0	1.5004×10−10	4.0
理论值		8		4

下载: 导出CSV 

| 显示表格

为对比本文算法与SPRD法^[11]的计算效果，对第三阶频率和振型的二次元解用SPRD法进行超收敛修复，用相邻四个结点(单元两端结点外加离单元最近的两个结点)上的位移解作三次多项式插值得单元上振型的超收敛解，将振型超收敛解代入Rayleigh商得频率超收敛解。频率和振型的收敛阶统计结果如表3所示。可见，SPRD法三次元修复解与本文算法三次元修复解的收敛阶相同，只是同样网格下SPRD法的解答精度比本文算法的解答精度稍低些。

表  3  例1第三阶特征对SPRD法收敛阶

Table  3.  Rate of convergence of recovered solutions on 3rd eigenpair in Example 1 with SPRD method

$ne$	$|\varOmega {\rm{ - }}{\varOmega ^{\rm{*}}}|$	$r$	$||{{d}} - {{{d}}^{\rm{*}}}||$	$r$
32	1.9168×10−5	−	6.1793×10−5	−
64	2.6890×10−7	6.2	3.9614×10−6	4.0
128	3.9339×10−9	6.1	2.4915×10−7	4.0
256	5.9333×10−11	6.1	1.5596×10−8	4.0
512	9.1038×10−13	6.0	9.7516×10−10	4.0
1024	1.4094×10−14	6.0	6.0954×10−11	4.0

下载: 导出CSV 

| 显示表格

例2. 常截面圆弧曲梁2

本例仍取轴线为半个圆周的常截面圆弧曲梁，曲梁参数同例1，两端支座形式改为面外简支，如图7所示。

图  7  例2常截面圆弧曲梁

Figure  7.  Uniform arch beam in Example 2

下载: 全尺寸图片 幻灯片

其振型可设为：

式中，$n \geqslant 0$且为整数。

将曲梁均匀划分为32个单元$(ne = {\rm{32)}}$，先用一次元$(m = {\rm{1)}}$进行有限元求解，再用二次元$(\bar m = {\rm{2)}}$进行超收敛修复。前12阶频率和振型的有限元解与超收敛解结果如表4所示。易见，对一次元，虽然振型结点位移不具有超收敛性，本文方法仍能显著提高解答(尤其是频率)的精度。

为考察本文算法对高阶频率和振型的计算效果，对该曲梁第26阶频率和振型进行分析，此阶对应$n = {\rm{14}}$。按前述方法计算该阶频率和振型的精确解，精确解为：

先用三次元$(m = {\rm{3)}}$进行有限元求解，再依次用四次元$(\bar m = 4)$、五次元$(\bar m = 5)$和六次元$(\bar m = 6)$进行超收敛修复，结果如图8和表5所示。结果表明p型超收敛算法对曲梁面外自由振动的高阶频率和振型仍有显著的修复效果。

表  4  例2前12阶特征对有限元解与超收敛解误差

Table  4.  Errors of FE solution and recovered solutions on first 12 eigenpairs in Example 2

阶数	$\varOmega$	${\varOmega ^h}$	$|\varOmega - {\varOmega ^h}|$	$||{{d}} - {{{d}}^h}||$	${\varOmega ^{\rm{*}}}$	$|\varOmega {\rm{ - }}{\varOmega ^{\rm{*}}}|$	$\dfrac{{|\varOmega - {\varOmega ^{\rm{*}}}|}}{{|\varOmega - {\varOmega ^h}|}}/( \text{%})$	$||{{d}} - {{{d}}^{\rm{*}}}||$	$\dfrac{{||{{d}} - {{{d}}^{\rm{*}}}||}}{{||{{d}} - {{{d}}^h}||}}/( \text{%})$	n
1	0.000000000	0.278903246	2.79×10−1	1.21×10−3	0.003516615	3.52×10−3	1.26	8.75×10−6	0.72	1
2	2.096810333	2.379666689	2.83×10−1	4.89×10−3	2.097006671	1.96×10−4	0.07	2.11×10−4	4.32	2
3	6.316663854	6.817878189	5.01×10−1	1.11×10−2	6.317460674	7.97×10−4	0.16	9.77×10−4	8.76	3
4	12.156595065	13.009811658	8.53×10−1	1.95×10−2	12.159147906	2.55×10−3	0.30	2.62×10−3	13.41	4
5	17.335703469	17.337172055	1.47×10−3	1.20×10−3	17.335703469	0.00×100	0.00	7.65×10−6	0.64	1
6	19.209142151	20.568557976	1.36×100	3.09×10−2	19.216004983	6.86×10−3	0.50	5.86×10−3	18.93	5
7	22.125304755	22.143921945	1.86×10−2	4.81×10−3	22.125316688	1.19×10−5	0.06	8.89×10−5	1.85	2
8	28.486387871	28.560421560	7.40×10−2	1.09×10−2	28.486493560	1.06×10−4	0.14	4.06×10−4	3.74	3
9	27.163609617	29.207708099	2.04×100	4.41×10−2	27.179901316	1.63×10−2	0.80	1.15×10−2	26.16	6
10	35.578254236	35.767316820	1.89×10−1	1.90×10−2	35.578729685	4.75×10−4	0.25	1.09×10−3	5.76	4
11	35.780085885	38.709492937	2.93×100	5.95×10−2	35.815415938	3.53×10−2	1.21	2.16×10−2	36.35	7
12	43.038408009	43.423183936	3.85×10−1	3.01×10−2	43.039902934	1.49×10−3	0.39	2.41×10−3	8.00	5

下载: 导出CSV 

| 显示表格

表  5  例2第26阶特征对有限元解与超收敛解收敛阶

Table  5.  Rate of convergence of FE solution and recovered solutions on 26th eigenpair in Example 2

$m = 3$					$\bar m = {\rm{4}}$
$ne$	$|\varOmega {\rm{ - }}{\varOmega ^h}|$	$r$	$||{{d}} - {{{d}}^h}||$	$r$	$ne$	$|\varOmega {\rm{ - }}{\varOmega ^{\rm{*}}}|$	$r$	$||{{d}} - {{{d}}^{\rm{*}}}||$	$r$
16	3.0092×10−1	−	8.6231×10−1	−	16	1.1674×10−2	−	8.4742×10−1	−
32	5.8019×10−3	5.7	2.5207×10−3	8.4	32	4.4443×10−5	8.0	2.0138×10−4	12.0
64	9.5518×10−5	5.9	1.4991×10−4	4.1	64	1.7980×10−7	7.9	5.2408×10−6	5.3
128	1.5120×10−6	6.0	9.2470×10−6	4.0	128	7.0914×10−10	8.0	1.5277×10−7	5.1
256	2.3701×10−8	6.0	5.7601×10−7	4.0	256	2.7769×10−12	8.0	4.6824×10−9	5.0
512	3.7063×10−10	6.0	3.5970×10−8	4.0	512	1.0854×10−14	8.0	1.4560×10−10	5.0
理论值		6		4	理论值		8		5
$\bar m = {\rm{5}}$					$\bar m = {\rm{6}}$
$ne$	$|\varOmega {\rm{ - }}{\varOmega ^{\rm{*}}}|$	$r$	$||{{d}} - {{{d}}^{\rm{*}}}||$	$r$	$ne$	$|\varOmega {\rm{ - }}{\varOmega ^{\rm{*}}}|$	$r$	$||{{d}} - {{{d}}^{\rm{*}}}||$	$r$
16	2.2943×10−3	−	8.4742×10−1	−	16	2.1079×10−3	−	8.4742×10−1	−
32	4.2573×10−7	12.4	1.4481×10−4	12.5	32	2.1331×10−7	13.3	1.4365×10−4	12.5
64	2.5609×10−10	10.7	2.4027×10−6	5.9	64	4.1349×10−11	12.3	2.3581×10−6	5.9
128	2.2096×10−13	10.2	3.8469×10−8	6.0	128	9.4404×10−15	12.1	3.7307×10−8	6.0
256	2.0927×10−16	10.0	6.0524×10−10	6.0	256	2.2647×10−18	12.0	5.8473×10−10	6.0
512	2.0281×10−19	10.0	9.4732×10−12	6.0	512	5.5045×10−22	12.0	9.1435×10−12	6.0
理论值		10		6	理论值		12		6

下载: 导出CSV 

| 显示表格

图  8  例2第26阶特征对有限元解与超收敛解收敛阶

Figure  8.  Rate of convergence of FE solution and recovered solutions on 26th eigenpair in Example 2

下载: 全尺寸图片 幻灯片

例3. 一般轴线曲梁

文献[27-28]研究了图9所示的抛物线、三角正弦线和椭圆线三类轴线曲梁的面外自由振动问题，三类曲梁轴线方程分别如下：

图  9  三类轴线曲梁

Figure  9.  Three types of curved beams

下载: 全尺寸图片 幻灯片

式中：$\xi = x/L$，$L$为曲梁跨度。三类轴线曲梁均采用相同参数，截面形状为矩形，$L/H = 5$，$L/b =$$28.87$，$t/b = 3$，$\nu = 0.3$，$k = 5/6$，$\varepsilon = 0.5$。记无量纲频率$\varOmega = \omega {L^2}\cdot$$\sqrt {\rho A/EI}$。

将三类轴线曲梁沿轴线水平投影方向均分为32个单元$(ne = 32{\rm{)}}$，先采用一次元$(m = 1)$进行有限元计算，后采用二次元$(\bar m = 2)$进行超收敛修复和先采用二次元$(m = 2)$进行有限元计算，后采用三次元$(\bar m = 3)$进行超收敛修复，计算前3阶频率如表6所示。与文献结果相比，频率在一次元有限元计算后仍误差较大，但采用二次元进行超收敛修复后，误差显著降低。而用二次元有限元计算后采用三次元进行超收敛修复，得到的解答与文献结果已十分接近。表7给出了上述网格下三类轴线曲梁前20阶振型有限元计算和超收敛计算的总耗时。可见，p型超收敛计算耗时很少，但解答精度提升非常显著。在一次元有限元计算后用二次元进行超收敛修复，所得结果与直接用二次元进行有限元计算的精度已很相近，但计算总耗时不足后者的1/2，充分展现出p型超收敛算法提升精度的高效性。

表  6  例3三类轴线曲梁的无量纲频率

Table  6.  Dimensionless frequencies for three curved beams of Example 3

轴线与支座	阶数	文献[27]	文献[28]	$m = {\rm{1}}$，$\bar m = {\rm{2}}$					$m = {\rm{2}}$，$\bar m = {\rm{3}}$
				有限元解	误差/(%)	超收敛解	误差/(%)		有限元解	误差/(%)	超收敛解	误差/(%)
抛物线铰支-铰支	1	6.0825	6.090	7.813	28.45	6.095	0.21		6.083	0.00	6.082	0.00
	2	30.4026	30.40	35.739	17.55	30.426	0.08		30.405	0.01	30.402	0.00
	3	70.0449	70.03	81.753	16.71	70.122	0.11		70.048	0.00	70.032	−0.02
三角正弦线固支-固支	1	17.0609	17.14	20.088	17.74	17.068	0.04		17.064	0.02	17.061	0.00
	2	48.5761	48.95	56.542	16.40	48.609	0.07		48.593	0.03	48.575	0.00
	3	95.0274	96.05	107.806	13.45	101.690	7.01		95.075	0.05	95.020	−0.01
椭圆线铰支-固支	1	11.0158	11.04	13.299	20.73	11.034	0.16		11.017	0.01	11.016	0.00
	2	38.6820	38.82	45.353	17.24	38.722	0.10		38.691	0.02	38.682	0.00
	3	81.9146	82.33	95.626	16.74	82.051	0.17		81.947	0.04	81.911	0.00

下载: 导出CSV 

| 显示表格

表  7  例3前20阶振型计算总耗时

Table  7.  Total computation time for first 20 modes of Example 3 /s

轴线	有限元解1 (m=1)	超收敛解1 ($\bar m = 2$)	有限元解2 (m=2)	超收敛解2 ($\bar m = 3$)
a)	12.204	0.063	30.223	0.105
b)	13.180	0.067	37.097	0.162
c)	14.292	0.084	38.639	0.146

下载: 导出CSV 

| 显示表格

综合数值算例，本文超收敛解的收敛规律如表8所示。

表  8  超收敛解的收敛阶

Table  8.  Convergence order of recovered solutions

收敛阶	${\omega ^*}$	${{{d}}^{\rm{*}}}$
$r$	$\min (2\bar m,4m)$	$\min (\bar m + 1,2m)$

下载: 导出CSV 

| 显示表格

5   结论

本文的p型超收敛算法在有限元解的基础上，充分利用频率和振型结点位移所具有的超收敛特性，通过逐单元求解振型近似满足的线性边值问题获得频率和振型的更高精度解，方法思路巧妙，稳定高效，能显著提高解答精度和收敛阶，可进一步推广应用到包括薄壁结构等空间结构自由振动的分析中。