孙浩涵，袁 驷

(1. 中物院高性能数值模拟软件中心，北京 100088；2. 北京应用物理与计算数学研究所，北京 100088； 3. 清华大学 土木工程系，北京 100084)

自由振动反映结构动力特性，是冲击等复杂力学分析的基础。在数学模型上，结构的自由振动可归结为偏微分方程特征值问题，通常以有限元法等数值方法进行求解，在各类实际工程分析中得到了广泛应用[1-2]。有限元法通用灵活，但作为一种近似数值方法会引入离散误差，网格划分将决定求解的效率和精度。自适应有限元法在传统有限元法(finite element method,FEM)的基础上，引入误差估计算法和网格细化技术，自适应地反复迭代和调整网格，最终提供优化的有限元网格和满足精度要求的解答。这一概念最早于20世纪70年代末由Babuška等[3-4]提出，首先被应用于线性问题，而后被推广至各类非线性问题。

特征值问题作为一类特殊的非线性问题，通常需要迭代求解，合理的网格划分将显著减少计算量。误差估计是特征值问题自适应分析的核心环节。在理论方面，Strang等[5]、Babuška等[6]、Knyazev等[7]先后给出了特征值问题的先验误差估计；在应用方面，更有价值的是后验误差估计，相关研究包括利用残差基于对偶性构造误差估计[8]、基于余能原理分析特征值上下界[9]、基于p型离散和Rayleigh商的Taylor展开形式构造误差估计[10]以及基于对应线性边值问题的误差估计[11-12]等。基于振型超收敛拼片恢复解，王永亮[13]提出了中厚圆柱壳自适应有限元分析策略，成功将自适应策略运用至自由振动问题。

作为一类超收敛计算方法，相较于经典SPR法[14-15]，单元能量投影法(element energy projection，EEP)不依赖超收敛点[16]，且恢复精度更高，可实现有限元解的逐点误差估计。针对结构自由振动，Sun等[17]已建立了一套基于EEP法的自适应分析策略，对频率和模态同时进行误差控制，可连续求解若干阶结果，给出满足精度要求的频率及逐点满足误差限要求的模态[18-19]。

在实际工程应用中，也存在另一类需求，即只要求频率满足精度，而并不关心模态误差大小。本文提出了频率的后验误差估计算法，继而建立了整体频率误差和局部模态误差的转换关系，最终实现了以频率误差控制为目标的自适应有限元分析策略。本方法的有效性在二阶Sturm-Liouville问题(简称SL问题)及弹性薄膜自由振动问题上得到了应用验证。

1 背景知识

1.1 研究对象

为方便一般性的讨论，本文沿用数学表达，自由振动问题中的频率和模态分别对应于特征值问题的特征值(平方根)和特征函数，不再加以区分。问题的微分方程形式为

Lu=λru

(1)

以及相应给定的边界条件。式中：L为定义的微分算子；λ为特征值；u为特征函数；r为质量项。

式(1)对应的能量泛函为

(2)

式中：a(·,·)为能量内积；(·,·)为常规的双线性型。由式(2)的一阶变分，得弱形式方程为

δπ(u)=a(u,δu)-λ·(u,δu)=0

(3)

1.2 单元能量投影法

作为一种有效的有限元超收敛计算算法，EEP法是本文特征值问题误差估计的核心。其基本思想是通过假设有限元法的能量投影定理

a(u-uh,vh)=0, ∀vh∈Sh

(4)

在单个单元上近似成立，推导出具有超收敛性质的EEP解u*。然而，这一性质仅对一维问题成立，对于二维乃至三维问题，应采用逐维离散、逐维恢复的方法。如图1所示，二维问题有限元离散可视作两次一维有限元离散过程的叠加，进而其超收敛解可通过两次应用一维EEP技术接续计算获得。该方法的有效性已在一系列二维及三维问题的应用中得到了确认[20-22]。

图1 二维有限元网格逐维离散Fig.1 D-by-D discretization of 2D problems

2 总体策略

2.1 自适应分析目标

本文研究对象为特征值，自适应分析最终目标是：在精确解λk(k=1,2,…,n)未知的情况下，获得充分好的有限元网格πk，使得有限元解满足

(5)

式中，Tol为事先给定的误差限。由于精确解未知，实际采用以下误差控制准则

(6)

2.2 整体和局部误差

在自适应有限元分析中，既需要整体误差控制标准，以确定何时停机；又需要局部误差控制标准，以驱动网格迭代。在本文中，由于整体误差以特征值进行衡量，而局部误差仍需以特征函数进行评估，因此需要建立关系，将二者有机地联系起来。

采用h型网格加密方式的自适应求解，收敛过程往往经历“均匀加密”到“局部加密”两个阶段，自适应收敛率约为最佳收敛率：对于一维问题，特征值收敛率为2m，特征函数收敛率为m+1；对于二维问题，特征值收敛率为m，特征函数收敛率为(m+1)/2。

(7)

式中，n(n=1,2)为问题维度。若假设新网格下特征值误差刚好满足误差限，则由式(7)关系可得特征函数最大误差限为

(8)

误差大于该值的单元均应细分。在实际算法中，为保证鲁棒性，误差限总取为历史路径的最小值，即

(9)

2.3 整体算法

本文自适应有限元分析总体策略可总结为：

步骤1有限元解uh。在当前网格下(初始网格用户给定)，求解有限元解。

步骤2EEP解u*。依照Gauss积分点分布确定采样点位置，计算超收敛位移及其导数。

步骤3误差检验。计算超收敛特征值，进行特征值后验误差估计，若满足式(6)则停机；否则，依据式(9)确定局部特征函数误差限，检查所有单元，对未满足误差要求的单元进行标记。

步骤4单元细分。采用基于四叉树结构的单元细分方案，对标记单元进行细分，返回步骤1。

图2给出了本文特征值问题自适应分析流程。相较于原算法，主要区别以虚线流程框标注，即特征值误差检验和超限单元标记。下文中，分别以缩写APEF和APEV表示先后两种不同的自适应方法(adaptive procedure guided by errors of eigen function / value)，后者即为本文所提出的算法。

图2 特征值问题自适应分析流程Fig.2 Flow chart for adaptive analysis of eigenvalue problem

3 算法细节

3.1 有限元离散

采用有限元法对式(1)进行分析。在本文中，对于一维问题，采用m次多项式单元；对于二维问题，采用双m次多项式单元，并结合基于四叉树结构的单元细分法进行网格更新[23]，导出分层结构的网格(如图3所示)。

图3 二维问题分层网格划分Fig.3 Hierarchical mesh refinement for 2D problem

(10)

式中：Ni和Nj为m次多项式形函数；dij为相应的结点位移；单元检验函数vh采用相同的插值形式。

确定单元及网格形式后，有限元求解归结为求uh∈Sh满足虚功方程

a(uh,vh)=λ·(uh,vh), ∀vh∈Sh

(11)

经由标准的矩阵集成流程，无论一维或二维问题，原问题式(1)转化为广义矩阵特征值问题

KD=λMD

(12)

式中：D为振型向量；K和M分别为整体刚度矩阵和一致质量矩阵。

3.2 特征值计算

式(12)的广义矩阵特征值问题，基于计数法和移位的逆幂(子空间)迭代，采用“划界”和“求解”两阶段法求解，可保证不丢根、不落根，并已在APEF算法中得到了充分验证；额外的精度导护措施，包括“μ检验”和“步长放大技术”，在本文中也同样采纳以确保解答的鲁棒性。

特征值通过计算Rayleigh商获得，可有效加速收敛过程。式(1)的各阶次解将使得Rayleigh商，即式(13)取为驻值

(13)

(14)

3.3 特征值误差估计

对式(13)取一阶变分，有

(15)

取驻值条件，令一阶变分为零，即得到式(3)。取二阶变分，并利用式(3)，有

(16)

由级数展开，有

λ-λh=δλ+δ2λ+O(δ3λ)+…

(17)

由于δλ=0，故有

(18)

其中，δu=eh=u-uh，则有特征值的误差估计

(19)

易知

a(eh,eh)≤C1h2m, (eh,eh)≤C2h2m+2

(20)

式中：h为单元尺寸参数；C1和C2为问题相关常数。故有Rayleigh商的误差估计

|λ-λh|≤Ch2m

(21)

式中，C为问题相关常数。可知有限元解的Rayleigh商误差为2m阶。

3.4 超收敛特征值计算

由式(21)，采用Rayleigh商计算特征值，其收敛阶主要取决于位移和导数精度。若可提供收敛阶更高的位移和导数，则理论上可获得超收敛的特征值，也即

(22)

不失一般性地，本文以二阶SL问题及弹性薄膜自由振动为模型问题进行讨论，数学表达式分别为

(23a)

(23b)

式中：n为边界外法向；ΓD为Dirichlet边界；ΓN为Neumann边界。

在算法设计时，有两个方面需要着重说明。

3.4.1 单元积分策略

在先前基于EEP法的特征值问题自适应有限元分析研究中，采样点主要用于单元上特征函数的误差估计，因而按照图4中采样方式(1)均匀分布；不同于此，本文采样点的另一作用在于作为积分点使用，以准确地计算式(22)，因而对于二维问题，应按照图4中采样方式(2)，也即相应数量Gauss点的位置分布。此外，对于本文采用的m次单元，通过EEP法获得的超收敛解可能由m+2次的多项式表达。为保证积分准确性，保证在各方向上至少有m+2个采样点。

图4 单元采样点分布Fig.4 Sampling points distribution over element domain

3.4.2 超收敛解的导数计算

在式(22)中，能量内积a(·,·)表达式中的积分项内一般包含位移的导数项。在本文中，采用对超收敛的EEP位移直接求导的方式进行计算。

对于SL问题，超收敛位移按式(24)计算，可直接求导

(24)

对于二维Laplace算子的特征值问题，两个方向的超收敛位移导数需要区分计算。本文中超收敛位移的计算采用单元逐维精度恢复方案，表达式为

(25)

(26)

图5 二维问题超收敛位移的导数计算Fig.5 Calculation of superconvergent derivatives for 2D problems

(27)

则对插值多项式PN-1(ξ)求导可得到ξ方向导数。

4 数值算例

本文算法已编制成Fortran 90程序。本章通过2个二阶SL问题算例、2个弹性薄膜自由振动算例，验证特征值超收敛计算及其驱动的自适应有限元分析的可靠性和有效性。

4.1 常截面杆轴向自由振动

本问题对应于式(23a)，参数p=r=1，q=0。对于该题，精确的特征值及特征函数表达为

λk=kπ2，uk=sin(kπx)

(28)

分别取低阶和中高阶的若干特征值，采用单元二分加密的方式，统计有限元解和超收敛解的收敛阶。

如表1～表3所示，对于常系数问题，用有限元解uh计算和用EEP解u*计算特征值，后者收敛阶可提升4阶，即由2m阶提升为2(m+2)阶。

表1 常截面杆轴向自由振动第1阶特征值收敛阶结果(线性元)Tab.1 The convergence order of the 1st eigenvalue of the axial free vibration of a bar with constant cross section (Linear elements)

表2 常截面杆轴向自由振动第1阶特征值收敛阶结果(三次元)Tab.2 The convergence order of the 1st eigenvalue of the axial free vibration of a bar with constant cross section (Cubic elements)

表3 常截面杆轴向自由振动第5阶特征值收敛阶结果(二次元)Tab.3 The convergence order of the 5th eigenvalue of the axial free vibration of a bar with constant cross section (Quadratic elements)

4.2 局部震荡问题

考虑式(29)参数定义的SL问题

p(x)=e10/cosh(10x),q(x)=0,
r(x)=100π2e10cosh(10x)/sinh210,
a=0,b=1,u(a)=0,u(b)=0

(29)

本例有精确解

(30)

该问题特点是在求解域(a,b)右端点附近的各阶特征函数均呈现剧烈震荡，求解有一定的困难。

首先验证超收敛特征值的有效性。表4给出了该例第3阶特征值的收敛阶，对于该变系数问题，用EEP解求征值其精度可提升2阶，即由2m阶提升为2(m+1)阶。

表4 震荡问题第3阶特征值收敛阶结果(二次元)Tab.4 The convergence order of the 3rd eigenvalue of oscillation problem (Quadratic elements)

设误差限为Tol=10-5，分别采用线性元及二次元，以APEV进行自适应分析。图6～图9对比了自适应及均匀网格加密的收敛过程，并绘制了单元尺寸分布。可以看到，基于本文自适应方案，收敛效率均有一定程度的提高，单元分布反映了问题的震荡特性。

图6 震荡问题自适应及均匀加密过程(线性元)Fig.6 Adaptive and uniform mesh refinement for oscillation problem (Linear elements)

图7 震荡问题自适应网格单元尺寸分布(线性元)Fig.7 Element size distribution of adaptive mesh for oscillation problem (Linear elements)

为进一步说明APEV的特点，对本例设置更苛刻的误差限Tol=10-8，作为对照，同时采用APEF进行计算。

采用三次元，从4个均分单元的初始网格开始，对第3阶特征对进行自适应分析。图10给出了两种方法自适应过程的误差下降图。可以看到，由于同时控制特征值和特征函数误差，APEF自适应方法存在较大的计算冗余度，算法停机时自由度数为947，特征值误差为1.01×10-11，远小于给定误差限；而本文算法直接控制特征值误差，算法停机时自由度数为308，特征值误差为7.06×10-9，体现了较高的计算效率。

图8 震荡问题自适应及均匀加密过程(二次元)Fig.8 Adaptive and uniform mesh refinement for oscillation problem (Quadratic elements)

图9 震荡问题自适应网格单元尺寸分布(二次元)Fig.9 Element size distribution of adaptive mesh for oscillation problem (Quadratic elements)

图10 震荡问题新旧方法自适应过程Fig.10 Adaptive mesh refinement with current and previous methods for oscillation problem

4.3 固支扇形膜自由振动

考虑四边固支的扇形膜(如图11所示)，其模态对应于环形弹性膜的反对称模态。本例中存在圆弧几何边界，具有一定的特殊性，采用精确几何单元捕捉该特点。

图11 扇形膜及三次几何精确单元Fig.11 Fan-shaped membrane and cubic geometrically exact element

本例存在精确解。对于半径为a≤r≤b的环形膜，其解析解为

(31)

式中，Jm和Ym分别对应于第一类和第二类的m阶Bessel函数。解析的特征值表达式为

(32)

式中，km,n为式(33)的第n阶根

(33)

对于本例，采用四次元，以单个单元网格为初始状态，进行前8阶特征值的自适应有限元分析，误差限设为Tol=10-4，结果汇总于表5。误差比显示，各阶特征值有限元解均满足误差要求。

表5 固支扇形膜自由振动特征值结果(APEV)Tab.5 Eigenvalues for free vibration of clamped fan-shaped membrane (APEV)

对本例，也采用APEF自适应策略进行分析，给定与以上相同的初始网格、单元及误差限。表6给出了计算结果。与APEV相比，APEF的冗余度较高，以第1阶为例，尽管使用了高达4 081个自由度，但其特征值结果反而差于APEV计算结果。

表6 固支扇形膜自由振动特征值结果(APEF)Tab.6 Eigenvalue for free vibration of clamped fan-shaped membrane (APEF)

图12和图13分别给出了第7阶的自适应最终网格及有限元解模态。

图12 固支扇形膜第7阶最终网格Fig.12 7th final mesh of clamped fan-shaped membrane

图13 固支扇形膜第7阶模态Fig.13 7th mode for clamped fan-shaped membrane

为更充分说明自适应网格划分的有效性，图14给出了均匀网格下第7阶模态的有限元解误差分布，与图12相一致地，自适应网格细密处误差也较大，体现了自适应单元加密的合理性。

以第4阶为例，图15对自适应过程中特征值的收敛性进行了分析。图15中，FEM和EEP分别表示在相同网格下经由有限元法和EEP后处理得到的特征值。在各个迭代步骤，基于超收敛特征值的误差估计均逼近真实误差，体现了其有效性；在停机时，误差小于误差限Tol，体现了APEV过程的有效性。

图14 固支扇形膜均匀网格下第7阶模态误差Fig.14 Errors of 7th mode for clamped fan-shaped membrane on uniform mesh

图15 固支扇形膜第4阶特征值收敛分析Fig.15 Convergence analysis of the 4th order eigenvalues of clamped fan-shaped membrane

4.4 固支L形膜自由振动

考虑固支L形膜自由振动问题(如图16所示)，其部分阶模态位移导数在凹角处存在奇异性，采用常规有限元法具有较高的求解难度。采用三次元，以三个单元的网格作为初始状态，误差限设为Tol=10-3。

图16 固支L形膜计算模型及初始网格Fig.16 Computational model and initial mesh for clamped L-shaped membrane

图17给出了第1、第3、第5阶的自适应最终网格及相应阶次模态，本文算法自动识别了相应阶次的模态特性。图18及图19给出了第1阶和第5阶特征值的收敛过程。对于本例，由于应力奇异点的存在，EEP法后处理特征值精度在自适应前期并不优于有限元解，但局部误差估计仍较为有效地体现了误差分布特性，体现在网格可自动识别问题难点，向奇异点附近加密，并以与误差限Tol相近的真实误差停机。

图17 固支L形膜计算模型及初始网格Fig.17 Meshes and modes for L-shaped membrane

图18 固支L形膜第1阶特征值收敛分析Fig.18 Convergence analysis of the 1st order eigenvalues of clamped L-shaped membrane

图19 固支L形膜第5阶特征值收敛分析Fig.19 Convergence analysis of the 5th order eigenvalues of clamped L-shaped membrane

5 结 论

本文建立了以频率(特征值)误差控制为目标的自适应有限元分析策略。以特征值误差估计为核心，本文首先分析了Rayleigh商的收敛阶，基于此提出了超收敛特征值的计算方案，其与原有限元解之差即可作为特征值的误差估计。以特征值误差估计作为整体停机标准，以特征函数误差估计衡量局部单元误差，本文通过两次网格迭代间的误差关系，将整体和局部的关系联系起来，最终形成了完整自适应方案。

自适应分析的主要目标是采用“尽可能少”的自由度数获得“尽可能高”的计算精度，其实现依赖于网格单元的合理分布，进一步依赖于对有限元解的有效误差估计。本文数值结果显示，对于光滑问题，相较于有限元解，EEP后处理特征值具有超收敛性，精度至少高2阶；对于奇异问题，特征值尽管不具备超收敛性，APEV自适应分析过程仍具有有效性，显著提升了有限元分析效率。此外，本文方法并不局限于作为模型的SL问题和弹性薄膜自由振动问题，在满足分层结构网格的条件下，该方法可一般性地推广至包括板、壳、三维弹性体在内的各类结构自由振动乃至弹性稳定问题。