平面曲梁面外自由振动有限元分析的p 型超收敛算法

2020-10-29叶康生

工程力学 2020年10期

叶康生，梁童

(清华大学土木工程系，土木工程安全与耐久教育部重点实验室，北京 100084)

平面曲梁结构在现代土木、航天、机械等领域具有广泛的应用，其自由振动是结构分析的重要内容。该问题目前的分析方法主要有：有限元法[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 所示，建立曲梁分析的xyz坐标系，x、y分别沿轴线切向、法向，z垂直于轴线所在平面。记s为轴线弧长坐标。曲梁面外自由振动时，分别记面外位移、y轴转角位移和扭转角位移为w(s) 、 θ(s) 、 ϕ(s) ，对应内力Q、M、T如图1 所示。

图1 平面曲梁Fig.1 Planar curved beam

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

式中： ()′=d()/ds；E为弹性模量；G为剪切模量； ρ为材料密度；A为截面面积；Iy为截面关于y轴的惯性矩；Ip为截面极惯性矩；k为截面形状系数；J为截面扭转常数；R为轴线曲率半径，当轴线曲率中心位于y轴正向时R为正，位于y轴负向时R为负。

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

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

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

记网格尺寸：

图2 有限元网格划分Fig.2 Finite element mesh

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

式中： ξ=(s−se)/he，为单元 (e)上的无量纲坐标；we、 θe、 ϕe为第e个结点上的结点位移；ai、bi、ci(i=1,2,···,m−1) 分别为w、 θ 、 ϕ的单元内部自由度。将式(9)写成矩阵形式：

式中， ∆e为单元 (e)的自由度向量，具体为：

Nw、Nθ、Nϕ分别为w、 θ 、 ϕ的形函数矩阵，其具体表达式可由式(9)和式(11)导出。

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

式中：Kh和Mh分别为该网格下的整体刚度矩阵和整体质量矩阵，分别由单元刚度矩阵和单元质量矩阵集成获得； ∆h为该网格下的结构整体自由度向量。单元刚度矩阵和单元质量矩阵具体计算如下：

求解上述矩阵广义特征值问题即可获得该网格下的有限元解 (λhn,dnh(s))，n=1,2,···。

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

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

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

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

2 超收敛求解思路

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

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

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

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

3 误差估计

4 数值算例

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

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

例1. 常截面圆弧曲梁1

图3 例1 常截面圆弧曲梁Fig.3 Uniform arch beam in Example 1

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

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

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

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

图4 例1 第三阶频率、振型和结点位移的收敛阶Fig.4 Rate of convergence of 3rd frequence, mode and its nodal displacement in Example 1

为考察p型超收敛算法的计算精度和修复效果，将曲梁均匀划分为5 个单元，先用一次元(m=1) 进行有限元求解，再用二次元 (=2)进行超收敛修复。得到第三阶频率的有限元解为Ωh3=8.398 (相对误差为300.5%)，超收敛解为Ω∗3=2.664(相对误差为27.0%)。图5 为第三阶振型的有限元解、超收敛解和精确解的对比。可见，本文算法可以显著提升频率和振型的计算精度。

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

解的收敛阶与=2m时相同，不会随的进一步提高而增加。

表1 例1 第三阶频率、振型和结点位移的收敛阶Table 1 Rate of convergence of 3rd frequency, mode and its nodal displacement in Example 1

图5 例1 第三阶振型的有限元解与超收敛解比较Fig.5 Comparison of FE solution with recovered solution on 3rd mode in Example 1

图6 例1 第三阶特征对有限元解与超收敛解收敛阶Fig.6 Rate of convergence of FE solution and recovered solutions on 3rd eigenpair in Example 1

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

例2. 常截面圆弧曲梁2

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

其振型可设为：

表2 例1 第三阶特征对有限元解与超收敛解收敛阶Table 2 Rate of convergence of FE solution and recovered solutions on 3rd eigenpair in Example 1

式中，n≥0且为整数。

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

表3 例1 第三阶特征对SPRD 法收敛阶Table 3 Rate of convergence of recovered solutions on 3rd eigenpair in Example 1 with SPRD method

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

图7 例2 常截面圆弧曲梁Fig.7 Uniform arch beam in Example 2

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

表4 例2 前12 阶特征对有限元解与超收敛解误差Table 4 Errors of FE solution and recovered solutions on first 12 eigenpairs in Example 2

例3. 一般轴线曲梁

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

图8 例2 第26 阶特征对有限元解与超收敛解收敛阶Fig.8 Rate of convergence of FE solution and recovered solutions on 26th eigenpair in Example 2

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

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