许多工程实际问题其求解区域往往要假设成无穷域，如土-结构相互作用、路面结构分析、渗流、波的传播、沉积物迁移等[1]。求解这些无穷域问题一般有两种方法[2]：一种是截取较大的有限区域进行计算，其存在两个缺陷，一是求解的计算量过大，二是无法满足实际问题在无穷远处的边界条件；另一种是在数值算法中使用无穷单元。
针对有限元法(Finite Element Method，简称FEM)这一通用的数值计算方法，无穷单元的研究和应用已取得一定的进展。1973年，Ungless[3]首次提出了无穷单元的概念，此后Bettess[4]提出乘子型无穷元，将有限元(Finite Element，简称FE)的形函数乘以一个衰减函数以插值得到衰减的位移。1983年，Zienkiewicz等[5]和Curnier[6]提出映射型无穷元，利用坐标变换函数，将无穷区域映射到有限区域上。基于这两种方法，无穷元被应用于土-结构相互作用[7]、波动[8－9]、热动力[10]、路面结构[11－12]等各类问题上。对于乘子型无穷元，一种衰减函数即一种无穷元，往往只对某类问题适用；无穷单元和相邻的常规有限单元之间有时会不协调；需要在半无穷域上进行数值积分，计算较为复杂，因此其应用范围不及映射型无穷元广泛。而映射型无穷元除坐标映射外，其余计算均与传统FEM相同，因而得到更多的青睐，但是其计算精度和效率较大地依赖于参与坐标映射的最远端点即极点的选择，还需要针对极点位置选择进行相关研究[13]。
除FEM的无穷元外，Beer和Watson [14]提出了边界元法中的无穷元，由于其求解基础是需要一个基本解，对于很难有基本解的复杂问题，方法受到局限。而我国学者也针对一些原创的新型数值计算方法提出了相应的无穷元[15－17]，并将其应用于实践。本文即围绕其中的一种——有限元线法(Finite Element Method of Lines，简称FEMOL)的无穷元展开。
FEMOL是一种半解析、半离散的数值计算方法。其基本思想是用一组结线将二维或三维求解区域离散成若干单元，基于变分原理并利用分片插值技术，得到以结线位移函数为未知量的常微分方程组(Ordinary Differential Equations，简称ODEs)，再用已有的通用ODEs求解器，如COLSYS[18]对其进行求解。如此便将偏微分方程问题转化成ODEs，对其建立映射型无穷元，将无穷边界的方向取作结线方向，则此结线方向上的解是由ODEs求解器得到，如COLSYS通过基于Gauss配点法的自适应求解得到，从而使得其结果对极点位置的依赖性大幅降低。可以说，FEMOL对无穷域问题具有天然的适应性。
本文在已有的FEMOL及其映射型无穷单元理论的基础上更进一步，将基于单元能量投影(Element Energy Projection，简称EEP)法的自适应FEMOL应用于无穷域问题的求解，用户只需输入稀疏的初始网格(一般包含一、两个单元即可)和期望满足的误差限Tl，就可以得到算法经误差估计、网格细分后自动生成的较为优化的最终网格，即自适应网格，在该网格上，FEMOL解的误差按最大模度量满足给定误差限。其中，无穷方向取为结线方向，由于自适应求解具有估计误差并调整网格划分的功能，因而无论极点位置优劣，最终结果都满足给定误差限。
文中首先介绍了二维FEMOL的基本原理、实施方法、无穷单元的构建，然后给出基于EEP法的自适应FEMOL策略，最后将其应用于圆柱绕流的Poisson方程问题、带孔无穷大板单向拉伸的弹性力学平面问题、半空间体受圆形均布荷载的三维轴对称问题，数值结果表明，该法用于求解二维线性无穷域问题取得了优良的效果，不仅能够给出满足误差限的函数(位移)解，也能给出具有优良性态的导数(应力)解。
1 FEMOL
本节以弹性力学平面问题为例简要介绍FEMOL的计算原理和方法。
1.1 常规FEMOL单元
二维问题中取x方向为FEMOL离散方向、y为结线方向，即将求解区域沿x方向进行网格划分，其中一典型单元经整体坐标到局部坐标的几何映射，转换为标准区域上的单元，如图1所示，其中Li表示第i条结线，Sj表示第j个端边，几何映射关系为：
这里
为m次Lagrange插值形函数。
图1 FEMOL二次单元的几何映射
Fig.1 A quadratic FEMOL element mapping
1.2 无穷结线与无穷单元
取无穷域问题的某个无限长方向为结线方向，对于拥有无限长结线(无穷结线)Li的单元即无穷单元，首先进行无穷结线的几何映射，如图2所示。
图2 无穷结线的映射
Fig.2 Infinite nodal line mapping
FEMOL中选用式(2)这一简单的映射关系：
式中：
为无穷结线的近端结点；
为无穷结线上的一个内部结点，
式(2)对应于图2映射的数学表达。一般称
为坐标变换的极点，并记a为近端结点到极点的距离。
无穷结线映射为有限结线之后，无穷单元可类似于图1的一般单元，几何映射到一个标准区域上，如图3所示。
图3 FEMOL二次无穷单元映射
Fig.3 A quadratic infinite FEMOL element mapping
这样，FEMOL中无穷单元和一般单元一样，转化为标准单元，后面的所有计算和处理即对标准单元实施开展。
1.3 FEMOL解
实施了上述线性坐标变换后，一般单元和无穷单元均转换为图1和图3中的标准单元，记为e。下面以弹性力学平面问题为例，给出FEMOL解的原理和算法。在其局部坐标系下，将e内任一点的位移u表示为对结线位移函数向量d e插值的形式：
式中，待求的单元结线位移向量d e、m次插值形函数矩阵N的表达如下：
式中：m为FEMOL单元位移插值次数，简称单元次数；
为Lagrange插值形函数；I2×2为2×2阶单位矩阵。
弹性力学平面问题的能量泛函为：
式中，
为作用在给定力边界Sσ上的分布力，能量内积a(u,v)定义为：
将式(2)代入式(6)，并将积分区间取为单元e，即可得到单元能量泛函Πe(de)，集成得到整体能量泛函
根据变分原理，由
可得到一组关于整体结线位移向量d的二阶ODEs：
以及相应的边界条件[19]。其中，系数矩阵A、G、H和向量F的所有元素均为已知量[19]。使用COLSYS求解式(8)，即在COLSYS自适应求解功能下得到满足其误差控制准则的结线位移解向量d，在各单元内由式(3)即可求得任意一点处的FEMOL位移解，记为uh。
2 基于EEP法的自适应求解
2.1 FEMOL的EEP超收敛算法
与FEM类似，由于采用分片插值技术，式(8)求得的结线位移解d精度较高，一般具有h2m的收敛阶，而由式(3)插值得到的内点位移解及其对ξ求导得到的ξ向导数，其精度较低，一般分别具有
的收敛阶[20]，且相邻单元间导数不连续。为了提高这些内点解的精度，并指导自适应分析，超收敛计算的研究由此展开[21―24]。
设检验函数为vh，精确解u和FEMOL解uh的误差记为
在整体求解域上易得精确满足的投影定理：
式中：
为检验空间。假设在单个FEMOL单元上，上述投影定理近似成立，即：
式中：
为线性形函数矩阵，
为线性形函数；
其中
式(9)这一假设成立的定理即称为单元能量投影(EEP)定理。由此定理，将单元e内局部坐标为
的一条直线记作
如图4所示，对于弹性力学平面问题，La位移的EEP超收敛计算公式[25]为：
当b[25]可逆时，单元内部任一点的EEP超收敛应力为：
当b不可逆时，单元内部任一点的EEP超收敛应力为：
式中：下标La表示线
上点的值；上标*表示EEP超收敛解；上标h表示FEMOL解；
表示u对应的斜面应力，
为几何映射的Jacobi值，b、c是与坐标和弹性常数有关的矩阵；下标ξ、η分别表示对ξ、η求导。即FEMOL超收敛计算的EEP法可以按式(12)～式(14)计算出整个二维求解域内任一点的超收敛的位移解以及应力解。
图4 La(ξ=ξa)位置示意图
Fig.4 Location of La(ξ=ξa)
EEP超收敛解具有很多优良的特性：1) 精度高：位移和导数的收敛阶均比常规FEMOL解至少高一阶；2) 应力连续：斜面应力在FEMOL相邻单元之间连续；3) 计算简便：是一种后处理计算方法，只要在常规FEMOL解上叠加一个积分项即可。由此在同类算法中体现出较大优势。
2.2 自适应求解策略
继续以弹性力学平面问题为例，本文自适应求解的目标是：在精确解u未知的情况下，预先给定误差限Tl，寻找一个优化的FEMOL网格划分π，使得在该网格上求得的FEMOL解uh的每个分量均按最大模度量逐点满足给定误差限Tl，即满足：
由于精确解u未知，以精度高于常规FEMOL解的EEP超收敛解u*代替真解u进行误差估计，即在最终网格π上要求所有单元均满足：
以此为本文自适应求解的误差估计方式和控制准则。按式(16)进行误差检验后，对于不满足的单元使用均差法[25]进行细分；直至所有单元均满足式(16)为止。
由此，可将基于EEP法的FEMOL自适应算法概括为如下三步。
1) FEMOL解：在当前网格上进行常规的FEMOL计算得到FEMOL解uh，初始时可取仅包含一、两个单元的稀疏网格作为当前网格。
2) EEP超收敛解：逐个FEMOL单元用式(12)计算其EEP超收敛位移解u*。
3) 网格细分：逐个FEMOL单元上，检验式(16)是否满足；对于所有未满足的单元，用均差法插入新的结线，得到新网格，设此网格为当前网格，返回第1)步；若所有单元都满足式(16)，则自适应求解结束，当前网格即为最终优化网格，此时的FEMOL解满足给定误差限Tl。
总的来说，上述EEP自适应技术具有以下优势：1) 算法程序自动给出适应问题物理和力学特性的优化网格，该网格下的位移解(即未知函数解)满足给定的精度要求，且一般较少冗余；2) 以位移(函数)解为求解目标，最终给出的常规FEMOL位移解按最大模度量逐点满足误差限，进一步地可由式(13)和式(14)求得精度较高的超收敛应力解；3) 算法思路明晰、实施简单，其基本思想也适用于FEM的自适应求解。
3 典型无穷域问题的自适应分析
3.1 自适应FEMOL对无穷域问题的适用性
由上述论述可见，常规单元和无穷单元经几何映射为标准单元后，FEMOL解、EEP超收敛解以及误差估计和网格细分都是直接对标准单元实施，因此研究上述自适应FEMOL对无穷域问题的适用性，关键在于其形如式(6)的能量泛函在无穷域上是否存在。
对于位移、应力均衰减的无穷域问题，该能量泛函一定存在，即为有限值，其中的能量内积a(u,v)也一定存在，这时无穷端边界直接按Dirichlet边界条件处理，且相应的投影定理式(9)也成立，单元投影定理式(10)和式(11)在无穷单元映射得到的标准单元上仍可假设为近似成立，EEP超收敛式(10)～式(14)适用，自适应FEMOL适用。对于能量泛函积分不收敛的情况，线性问题范畴内可通过分解叠加等策略对原问题稍作处理，对其中泛函收敛的部分应用自适应FEMOL求解，下文将针对具体问题加以说明。
在已有FEMOL常规单元自适应分析程序的基础上，本文以Fortran90语言编制了常规/无穷单元自适应FEMOL程序。对于具有代表性的二维无穷域问题，下文将给出计算结果，进一步阐述算法实施技术，并展示其优良效果。
3.2 问题1：带圆孔无限大板的单向拉伸
开孔构件在土木，机械，航空航天等工程领域有着广泛应用。很多情况下孔口尺寸会比构件尺寸小很多，直接用有限区域计算所需计算量较大，常常将其视为无穷域问题来分析。设无限大板上圆孔半径r0=0.1 m，单向均匀拉伸荷载q=20 GPa，板的弹性模量E=200 GPa，泊松比v=0.3。
显然，该问题的位移解与应力解都不是衰减的，因此直接求解时其能量泛函不收敛。采用分解叠加的策略可以解决此问题。如图5所示，将有孔受外荷载的原问题解分解为状态(I)：无孔单向拉伸板的解与状态(II)：有孔板在孔边作用有状态(I)孔边应力反向力的解两部分，状态(I)与状态(II)分别求解后再叠加。
对于状态(I)的求解十分容易，以孔心对应位置为坐标原点，在直角坐标下σx=q，σy=0，τxy=0。进而可得位移解为：
图5 带孔无限大板的叠加法求解示意图
Fig.5 Superposition for infinite plate with a circle hole
对于状态(II)，位移和应力都是衰减的，可以采用本文方法进行FEMOL自适应分析。由对称性，取1/4结构计算，FEMOL初始网格如图6所示，共2个单元。
图6 问题1的FEMOL初始网格
Fig.6 Initial FEMOL mesh of problem 1
首先在此初始网格上，无穷结线映射时取不同的极点位置进行常规FEMOL解，将计算结果与理论解[26]对比可以看出，极点位置a的取值对求解精度的影响很小，如表1所示，优于一般常规FEM[13]。
表1 不同极点位置下常规FEMOL解与理论解的误差对比
Table 1 Comparisons of error of FEMOL and exact solution among different pole points
a σ/(x=1, y=0) 误差/(%)/yq σ(x=2, y=0) 误差/(%)yq 0.1 2.9383 2.06 1.1742 3.66 0.5 3.0002 0.00 1.1775 3.38 1 3.0084 0.28 1.1784 3.31 10 3.0150 0.50 1.1817 3.04 100 3.0196 0.65 1.1787 3.31理论解 3.0000 ― 1.21875 ―
下面进行自适应求解。误差限给定为Tl=10-7 m。采用m=3次元，最终得到如图7所示的15个单元的FEMOL网格。
图7 问题1自适应求解的最终网格
Fig.7 The final adaptive mesh of problem 1
将最终自适应网格下的解叠加上状态(I)的解就得到了原问题的解。根据极坐标系(r,θ)下原问题的应力理论解[26]给出，由此推出原问题的位移解为：
将本文解与式(18)～式(19)的理论解相比较，其误差分布图如图8、图9所示，由于无法画出无限大的区域，这里截取2 m×2 m的范围，相对于圆孔半径0.1 m，足够可以展示本文解答的可靠性。显而易见域内任一点的两个位移分量的解答，其误差按最大模度量都满足给定误差限的要求。
进一步地，在自适应最终网格上由式(13)和式(14)求得超收敛应力解，叠加状态(I)的解答后，将其与理论解进行对比，图10和图11分别给出了r=0.1 m、θ∈(0,π/2)以及θ=π/4、r∈(0.1,1)m这两个区域内的EEP超收敛应力解与理论解的对比。从结果可以看出，尽管本文的自适应求解是以位移为目标的，但是最终网格上的EEP超收敛应力解不仅精度高，而且在FEMOL相邻单元之间连续，在应力边界处(r=0.1m)自动满足平衡条件，这是常规FEM或FEMOL无法做到的。
图8 问题1自适应FEMOL解uh(左)的误差分布图
Fig.8 The distribution of errors of adaptive FEMOL solutions uh (left) for problem 1
图9 问题1自适应FEMOL解与vh(右)的误差分布图
Fig.9 The distribution of errors of adaptive FEMOL solutions vh (right) for problem 1
图10 在r=0.1 m处EEP超收敛应力与理论解的比较
Fig.10 Comparisons of EEP stresses with exact stresses at r=0.1 m
图11 在θ=π/4处EEP超收敛应力与理论解的比较
Fig.11 Comparisons of EEP stresses with exact stresses at θ=π/4
3.3 问题2：半空间体作用圆形均布荷载
作用有圆形均布荷载的半空间体问题是典型的三维轴对称问题，可取如图12所示的任一对称面作为求解区域，将弹性力学平面问题中的坐标系(x,y)换成坐标系(r,z)进行计算，其中r为径向坐标，z为半空间体深度方向的坐标。
图12 问题2的求解区域示意图
Fig.12 Solution domain of problem 2
设半空间体材料的弹性模量E=100 MPa，泊松比取为0.3；圆形荷载作用半径为0.1 m，大小为1 MPa。FEMOL初始网格如图13所示，共3个单元，其中包含两个三角形无穷单元。可利用集中荷载作用在半空间体下的解析解[26]进行数值积分求得该问题的理论解。
图13 问题2的自适应FEMOL初始网格
Fig.13 Initial FEMOL mesh of problem 2
自适应求解的误差限给定为Tl =10-6m，采用m=3次元，最终得到如图14所示的14个单元的自适应网格。其最小单元的宽度为0.0008 m，紧邻荷载作用外边缘r=0.1 m处。越邻近此边缘，网格划分越密集，方法自动适应了问题的力学性质。
图14 问题2的自适应FEMOL求解最终网格
Fig.14 The final adaptive FEMOL mesh of problem 2
自适应网格下位移的FEMOL解和理论解相比的误差分布如图15、图16所示。由于无法画出无限大的区域，这里截取2 m×2 m的范围，相对于荷载作用范围0.1 m，足够可以展示本文解答的可靠性。可见，自适应网格下FEMOL位移解逐点具有满足误差限Tl的精度。
虽然自适应计算时只对位移进行了控制，但在自适应网格下计算出的超收敛应力解也有很高的精度。图17与图18给出了r=0 m和z=0 m处自适应网格下EEP超收敛应力解与理论解的比较。
图15 问题2中自适应FEMOL解vh的误差分布图
Fig.15 The distribution of errors of vh for problem 2
图16 问题2中自适应FEMOL解wh的误差分布图
Fig.16 The distributions of error of wh for problem 2
图17 在r=0 m处EEP超收敛应力解与理论解的比较
Fig.17 Comparisons of EEP stresses with exact stresses at r=0 m
图18 在z=0 m处EEP超收敛应力解与理论解的比较
Fig.18 Comparisons of EEP stresses with exact stresses at z=0 m
3.4 问题3：圆柱绕流
圆柱绕流问题是流体力学中的经典问题，指的是流体绕不转动的圆柱体的流动，如图19所示。以速度势φ为基本未知量，其控制方程与边界条件为：
图19 圆柱绕流示意图
Fig.19 Flow around a circular cylinder
由于常速度场的原因，该问题的速度势与速度都不是衰减的。这里再次采用分解叠加的策略进行处理。设常速度场对应的速度势为φ1，则对应于速度势φ2=φ-φ1的解答就可以使用FEMOL无穷单元的自适应分析来得到。
易得对应于速度势φ1的解为φ1=v0rcosθ，则φ2对应的控制方程与边界条件为：
现来求解φ2，由对称性，可取1/4结构计算。设速度v0=1 (m/s)、圆柱半径r0=1 m，FEMOL初始网格如图20所示，共2个单元。
自适应求解的误差限给定为Tl =10-4m ，采用m=3次元，最终得到如图21所示的10个单元的自适应网格。
图20 问题3 FEMOL初始网格
Fig.20 Initial FEMOL mesh of problem 3
图21 问题3最终自适应网格
Fig.21 The final adaptive mesh of problem 3
将最终自适应网格下的解叠加上常速度场的解就得到了原问题的解。原问题的理论解为[27]：
图22给出了速度势φ自适应FEMOL的解与理论解的误差分布，结果再次表明，求解域内任一点的FEMOL解的误差按最大模度量满足给定的误差限，本文方法可靠有效。
图22 问题3自适应FEMOL解φh的误差分布图
Fig.22 The distribution of errors of φh for problem 3
4 结论
FEMOL映射型无穷单元在求解无穷域问题时，其精度对极点位置的依赖性较小，有着很好的适用性。本文在已有的FEMOL及其映射型无穷单元理论的基础上，将具有无穷元的EEP法自适应FEMOL求解策略应用于典型的二维无穷域问题求解。数值结果表明，该法高效可靠，自适应分析能够给出满足误差限的函数(位移)解，在自适应网格上的超收敛导数(应力)解也性态优良。该法在三维问题、非线性问题等更广阔领域的应用值得进一步研究。
[1]Erkal A, Laefer D F, Tezcan S.Advantages of infinite elements over prespecified boundary conditions in unbounded problems [J].Journal of Computing in Civil Engineering, 2014, 29(6): 04014085.
[2]赵崇斌, 张楚汉, 张光斗.用无穷元模拟拱坝地基[J].水利学报, 1987, 18(2): 21－36.Zhao Chongbin, Zhang Chuhan, Zhang Guangdou.Simulation of arch dam foundation by using infinite element [J].Journal of Hydraulic Engineering, 1987,18(2): 21－36.(in Chinese)
[3]Ungless R F.Infinite finite element [D].Vancouver:University of British Columbia, 1973.
[4]Bettess P.Infinite elements [J].International Journal for Numerical Methods in Engineering, 1977, 11(1): 53－64.
[5]Zienkiewicz O C, Emson C, Bettess P.A novel boundary infinite element [J].International Journal for Numerical Methods in Engineering, 1983, 19(3): 393－404.
[6]Curnier A.A static infinite element [J].International Journal for Numerical Methods in Engineering, 1983,19(10): 1479－1488.
[7]Medina F.Modelling of soil-structure interaction by finite and infinite elements [J].Dissertation Abstracts International, 1980, 42(7), Section: B: 2929.
[8]Honjo Y, Pokharel G.Parametric infinite element for seepage analysis [J].International Journal for Numerical and Analytical Methods in Geomechanics, 1993, 17(1):45－66.
[9]Yang Y B, Hung H H, Lin K C, et al.Dynamic response of elastic half-space with cavity subjected to P and SV waves by finite/infinite element approach [J].International Journal of Structural Stability and Dynamics, 2015, 15(7): 1540009.
[10]Kazakov K S.Elastodynamic infinite elements with united shape functions for soil-structure interaction [J].Finite Elements in Analysis and Design, 2010, 46(10):936－942.
[11]Wang Q, Brill D R.Improvements in the application ofinfinite elements to the NIKE3D program for airport pavement response analysis [J].International Journal of Pavement Engineering, 2013, 14(5): 429－439.
[12]Liu P, Wang D, Oeser M.Application of semi-analytical finite element method coupled with infinite element for analysis of asphalt pavement structural response [J].Journal of Traffic and Transportation Engineering(English Edition), 2015, 2(1): 48－58.
[13]Abdel-Fattah T T, Hodhod H A, Akl A Y.A novel formulation of infinite elements for static analysis [J].Computers & Structures, 2000, 77(4): 371－379.
[14]Beer G, Watson J O.Infinite boundary elements [J].International Journal for Numerical Methods in Engineering, 1989, 28(6): 1233－1247.
[15]石绍春.有限元线法无穷单元[D].北京: 清华大学,1992.Shi Shaochun.Infinite FEMOL elements-Formulation and applications [D].Beijing: Tsinghua University, 1992.(in Chinese)
[16]袁驷, 石绍春, 崔京浩.有限元线法的无穷单元─无穷线的映射[J].数值计算与计算机应用, 1998,19(2): 81－89.Yuan Si, Shi Shaochun, Cui Jinghao.An infinite FEMOL element-Infinite line mapping technique [J].Journal on Numerical Methods and Computer Applications, 1998,19(2): 81－89.(in Chinese)
[17]袁帅, 钟宏志.无穷域问题的弱形式求积元分析[J].岩土力学, 2016, 37(4): 1187－1194.Yuan Shuai, Zhong Hongzhi.Analysis of unbounded domain problems by the weak form quadrature element method [J].Rock and Soil Mechanics, 2016, 37(4):1187－1194.(in Chinese)
[18]Ascher U, Christiansen J, Russell R D.Algorithm 569,COLSYS: Collocation software for boundary value ODEs [J].ACM Trans Math Software, 1981, 7(2): 223－229.
[19]Yuan S.The finite element method of lines: theory and applications [M].Beijing-New York: Science Press,1993.
[20]庞之垣.有限元线法的误差估计[J].计算数学, 1993,15(1): 110－120.Pang Zhiyuan.The error estimate of finite element method of lines [J].Mathematica Numerica Sinica, 1993,15(1): 110－120.(in Chinese)
[21]袁驷, 王枚, 王旭.二维有限元线法超收敛解答计算的EEP法[J].工程力学, 2007, 24(1): 1－10.Yuan Si, Wang Mei, Wang Xu.An element-energyprojection method for super-convergent solutions in two-dimensional finite element method of lines [J].Engineering Mechanics, 2007, 24(1): 1－10.(in Chinese)
[22]袁驷, 方楠, 王旭, 等.二维有限元线法自适应分析的若干新进展[J].工程力学, 2011, 28(3): 1－8.Yuan Si, Fang Nan, Wang Xu, et al.New progress in self-adaptive analysis of two-dimensional finite element method of lines [J].Engineering Mechanics, 2011, 28(3):1－8.(in Chinese)
[23]袁驷, 王永亮, 徐俊杰.二维自由振动的有限元线法自适应分析新进展[J].工程力学, 2014, 31(1): 15－22.Yuan Si , Wang Yongliang, Xu Junjie.New progress in self-adaptive FEMOL analysis of 2D free vibration problems [J].Engineering Mechanics, 2014, 31(1): 15－22.(in Chinese)
[24]徐俊杰, 袁驷.有限元线法EEP超收敛计算简约格式的再简约[J].工程力学, 2016, 33(增刊1): 1－5.Xu Junjie, Yuan Si.Further simplification of the simplified formulation of EEP super-convergent calculation in FEMOL [J].Engineering Mechanics, 2016,33(Suppl1): 1－5.(in Chinese)
[25]方楠.基于EEP法的弹性力学问题有限元线法自适应分析[D].北京: 清华大学, 2011.Fang Nan.Adaptive FEMOL analysis of elastic problems based on EEP super-convergent method [D].Beijing:Tsinghua University, 2011.(in Chinese)
[26]Timoshenko S P, Goodier J N.Theory of elasticity [M].New York: McGraw-Hill, 1987.
[27]郑邦民.流体力学[M].北京: 水利水电出版社, 1989.Zheng Bangmin.Hydromechanics [M].Beijing: China Water Power Press, 1989.(in Chinese)