变截面变曲率梁振型的有限元超收敛拼片恢复解和网格自适应分析

王永亮

(1.中国矿业大学(北京)力学与建筑工程学院，北京100083；2.中国矿业大学(北京)煤炭资源与安全开采国家重点实验室，北京100083)

梁构件在实际工程中广泛应用，如结构工程中的曲线梁桥[1]、航空航天工程中的机翼翼梁[2]、采矿工程中的矿柱Euler–Bernoulli梁模型[3]等。梁的自由振动是强迫振动研究的基础，自振频率可为避免结构共振提供依据[4]，利用振型叠加可分析地震载荷作用下结构的变形量[5]。此外，利用自振频率和振型可以进行含损伤结构的损伤识别[6−7]，损伤个数、位置和程度的准确识别依赖于高精度的频率和振型解答[8−9]。本文针对典型的平面变截面变曲率梁构件自由频率和振型求解开展研究，该类曲梁的轴线在同一平面内、且曲梁的任意位置横截面均关于上述平面对称。这类曲梁在面内振动和面外振动相互解耦，需要分别求解[10]，如何精确、有效地求解变截面变曲率梁的连续阶频率和振型成为需求，也是本文的研究目标。

有限元法是求解曲梁自由振动频率和振型近似解的重要方法，被应用于分析不同横截面[11]、支撑类型[12]、曲线线型形式[13]等复杂梁构件的自由振动问题。以控制和提高有限元解的精度为目标，自适应有限元法被应用于Euler-Bernoulli梁构件[6]、梁系结构[14−15]的频率和振型的求解，各类新型网格自适应算法[16− 17]为提供满足预设误差限的高精度解答提供了高性能计算方案。本文将建立变截面变曲率梁面内和面外自由振动振型的有限元超收敛拼片恢复解的计算方法，并基于振型超收敛解构建网格自适应分析策略，可获得优化的网格和满足预设误差限的连续阶频率和振型。数值算例检验了本文方法分析不同曲线形态、多类边界条件、变曲率、变化截面形式的曲梁面内和面外自由振动连续阶频率和振型的求解效力，表明本文方法分析结果的精确性和有效性。

1 曲梁面内和面外自由振动

考虑图1所示的平面变截面变曲率梁，曲梁的中性轴坐标为s，坐标系为xy z，其中x、y为曲梁平面内坐标，x沿轴线切向，y沿轴线法向，z垂直于轴线所在平面。面内振动的位移为：沿x轴位移振幅u、沿y轴位移振幅v和绕z轴的转角振幅ψz；面外振动的位移为：沿z轴位移振幅w、绕y轴的转角振幅 ψy和绕x轴的扭转角振幅φx。记曲梁曲率半径为R(s)，截面剪切刚度修正系数为κ，截面面积为A(s)，对y轴惯性矩为Iy(s)，对z轴惯性矩为I z(s)，极惯性矩为Ip(s)，扭转常数为J(s)，长度为l。记材料弹性模量为E，剪切模量为G，泊松比为ν，密度为 ρ。

图1 平面变截面变曲率梁坐标系和符号Fig.1 Coordinatesystemsand symbols of planar nonuniform and variablecurvature curved beam

本文研究的平面曲梁面内自由振动的微分控制方程[18]为：

式中：()′=d()/ds；ω为自振频率。

面内、面外自由振动控制方程式(1)、式(2)均可记为如下矩阵形式的特征值方程：

式中：u为对应的振型(位移)函数向量(面内、面外振动分别为{u,v,ψz}T、{w,ψy,φx}T)，ω、u分别对应于特征值、特征向量，本文亦将(ω,u)合称为特征对；L、R是相应的微分算子矩阵。

2 振型有限元解

对于求解特征值方程式(3)，在给定有限元网格π下，常规有限元建立如下线性矩阵特征值方程：

式中：D为振型向量的有限元解；K和M是静力刚度矩阵和一致质量矩阵。采用逆幂迭代法[20]求解，即得当前网格下有限元解(ωh,u h)。

3 振型超收敛拼片恢复解

位移型有限元法以节点位移为基本未知量，有限元位移解在单元内部具有hm+1阶(m为当前形函数阶次)的超收敛性，在单元节点上具有h2m阶的超收敛性，这些节点成为位移超收敛点。超收敛点上解答的误差比其他区域解答更小，具有更快的收敛速度[20]。利用相邻单元拼片，并将这些单元节点位移值进行次高阶形函数(高于当前形函数阶次m，即≥m+1)插值，即可获得单元域内位移场的超收敛解，形成有限元后处理超收敛拼片恢复方法[6,21−22]。在当前网格下，位移超收敛解比常规有限元解具有更高收敛阶。本文对于曲梁的面内和面外自由振动问题，求得当前网格下振型(位移)的有限元解后，利用有限元后处理超收敛拼片恢复方法，将单元e的相邻单元进行拼片，并进行高阶形函数插值处理，可以得到振型的超收敛解：

式中，a( )、b( )为应变能、动能内积。

运用Rayleigh 商计算得到的频率超收敛解比振型超收敛解具有更高的收敛阶[23]，因此，本文针对振型进行误差估计和控制，即可获得高精度的振型解，并确保频率解的准确性。

4 误差估计和网格细分加密

利用振型误差估计，网格可以进行优化处理来降低和控制振型的误差，达到预设的解答精度。本文方法采用式(14)对每个有限元单元e上的振型误差进行判断，如果误差控制式(14)不满足，则表明该单元上振型解答的误差过大，需要通过进行网格优化处理，本文采用单元均匀细分加密的方式来增加模型自由度、降低单元上解答的误差[6,24]。当前单元细分生成的新单元长度与目前误差和单元阶次相关，即利用当前误差可以估计新单元长度：

为适应问题的复杂性，往往在求解域上使用非均匀分布网格才能高效地获得可靠解答。上述自适应分析过程将单元均匀细分加密，如果每次细分出过多单元，则容易造成的最终单元的冗余，增大计算规模。本文方法对每次单元细分加密数目进行控制，单元均匀细分的实际数目为：

5 自适应分析策略

对曲梁自由振动问题，利用上述基于振型超收敛拼片恢复解的有限元解误差估计、网格细分加密方法，可以形成有效的网格自适应分析方法，最终可以得到优化的网格，并获得该优化网格上满足误差限的各阶高精度频率和振型解答。对于振型，超收敛解u∗较当前网格而下有限元解u h具有更高的收敛阶，即u∗比u h更接近精确解u，因此采用u∗代u对u h进行误差估计和控制，形成能量模形式的误差估计；对于频率，利用位移超收敛解计算Rayleigh 商[23]得到频率的超收敛解。通过Rayleigh 商得出的频率超收敛解比振型超收敛解具有更高的收敛阶，本文通过估计和控制振型解答的误差可保证频率解答满足误差要求，如此可形成本文算法通过控制振型误差的停机准则。基于有限元解的误差估计进行网格细分加密，即可求得优化的网格和满足误差限的各阶高精度频率和振型解答，形成一套有限元自适应求解方案：

1)有限元解：对于曲梁振动问题，在当前网格π 上进行常规有限元计算，得到当前网格下特征对的有限元解(ωh,u h)；

2)振型超收敛解：利用有限元后处理超收敛拼片恢复方法，可得当前网格下振型的超收敛解u∗；同时，可得超收敛解与有限元解之间的误差估计值ξ；

3)网格细分加密：对于 ξ不满足式误差控制式(14)的单元，用网格细分加密方法将其细分，得到新的有限元网格，返回步骤1)；如果所有单元均满足误差限，则网格无需再细分，求解过程结束。

6 数值算例

本文方法已经编制Fortran 90程序，本节给出求解多种变截面变曲率梁面内和面外自由振动的数值算例，通过对比典型的抛物线曲梁、变截面变曲率梁、椭圆弧曲梁面内自由振动和抛物线曲梁、圆弧曲梁面外自由振动结果来验证本文方法计算结果的精确性和适用性。本节所有算例均采用3次元，初始网格采用2 个单元，给定的初始误差限为Tol=10−4，设定单元细分的极限数目为d=6。

例1.抛物线曲梁面内自由振动

考虑如图2所示的抛物线曲梁，曲梁跨度L=5 m、高度h=4 m。

图2 抛物线曲梁几何模型Fig.2 Geometric model of parabolic curved beam

该曲线梁的抛物线方程为：

该曲线梁的基本参数如下：

表1 抛物线曲梁面内自由振动无量纲频率值Table 1 Non-dimensional frequencies of in-plane vibration of parabolic curved beam

图3 本文方法解答与密集网格高精度解的误差Fig.3 Errorsof results of proposed method and highprecision solutionsunder dense mesh

图4～图6分别给出了本文方法求解两端固支边界、高跨比h/L=0.8情况下的前3阶振型结果，并在横坐标轴上标记出自适应网格的最终分布情况。为方便直观显示和对比分析，图中振型结果均进行归一化处理(令最大振型值为1)。可以看出，振型在两固定端的位移均为0值；本文方法求解各阶振型均划分出非均匀网格，且在振型变化平缓区域使用稀疏网格、在振型变化剧烈处采用了相对细密的网格，避免了全域使用一致细密网格的冗余性。

图4 第1阶振型(uh1,vh1,ψ hz1)Fig.4 First order vibration mode(uh1,vh1,ψ hz1)

图5 第2阶振型(uh2,vh2,ψhz2)Fig.5 Second order vibration mode (uh2,vh2,ψ hz2)

图6 第3阶振型Fig.6 Third order vibration mode (uh3,vh3,ψ hz3)

例2.变截面变曲率梁面内自由振动

考虑图7所示的变截面变曲率梁，曲梁两端为固支边界条件、跨度为L。

图7 变截面变曲率梁几何模型Fig.7 Geometric model of non-uniform and variable curvaturecurved beam

使用本文方法连续求解该曲梁面内自由振动的连续前4阶特征对，计算频率值列于表2。文献[26]发展基于Wittrick-Williams算法的动力刚度法求解该问题，解答如表2所示，可以看出本文求解结果与该动力刚度法结果吻合较好，检验了本文方法对同时变化截面和曲率形式梁的面内自由振动求解的有效性。

表2 变截面变曲率梁面内自由振动频率值Table 2 Frequencies of in-plane vibration of non-uniform and variable curvature curved beam

例3.椭圆弧曲梁面内自由振动

例4.抛物线曲梁面外自由振动

考虑例1中图2所示抛物线曲梁的面外自由振动，抛物线方程同式(18)，该曲线梁的基本参数为：

图8 椭圆弧曲梁几何模型Fig.8 Geometric model of of elliptic curved beam

表3 椭圆弧曲梁面内自由振动无量纲频率值Table 3 Non-dimensional frequenciesof in-plane vibration of elliptic curved beam

表4 抛物线曲梁面外自由振动无量纲频率值Table 4 Non-dimensional frequencies of out-of-plane vibration of parabolic curved beam

表5 圆弧曲梁面外自由振动无量纲频率值Table5 Non-dimensional frequenciesof out-of-plane vibration of circular curved beam

例5.圆弧曲梁面外自由振动

考虑图10所示两端固支圆弧曲梁，圆弧角度为 θ0，梁横截面为圆形。

图10 圆弧曲梁几何模型Fig.10 Geometric model of circular curved beam

该曲线梁的基本参数为：

7 结论

本文研究工作的主要结论和计算方法潜力可总结如下：

(1)本文提出变截面变曲率梁面内和面外自由振动振型的有限元后处理超收敛拼片恢复方法，建立各阶振型的超收敛解，基于振型超收敛解进行网格自适应加密，求解得到优化的网格与满足用户给定误差限的频率和振型。

(2)本文算法具有通用性，分析了不同曲线形态、多类边界条件、变截面、变曲率形式的曲梁振动问题，可推广到一般形式结构和固体变形场有限元解的超收敛计算及自适应分析。

(3)本文方法具有应用分析含裂纹损伤曲线形杆件振动和稳定等工程实际问题的潜力，获得损伤缺陷下的高精度振型和失稳模态，为分析裂纹损伤对曲线形杆系结构的影响提供可靠解答。