留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于EEP技术的一维有限元结点位移误差计算

袁驷 邢沁妍 袁全

袁驷, 邢沁妍, 袁全. 基于EEP技术的一维有限元结点位移误差计算[J]. 工程力学, 2020, 37(9): 1-7, 29. doi: 10.6052/j.issn.1000-4750.2020.07.0446
引用本文: 袁驷, 邢沁妍, 袁全. 基于EEP技术的一维有限元结点位移误差计算[J]. 工程力学, 2020, 37(9): 1-7, 29. doi: 10.6052/j.issn.1000-4750.2020.07.0446
Si YUAN, Qin-yan XING, Quan YUAN. CALCULATION OF ERRORS OF NODAL DISPLACEMENTS IN ONE-DIMENSIONAL FINITE ELEMENT METHODS USING ELEMENT ENERGY PROJECTION TECHNIQUE[J]. Engineering Mechanics, 2020, 37(9): 1-7, 29. doi: 10.6052/j.issn.1000-4750.2020.07.0446
Citation: Si YUAN, Qin-yan XING, Quan YUAN. CALCULATION OF ERRORS OF NODAL DISPLACEMENTS IN ONE-DIMENSIONAL FINITE ELEMENT METHODS USING ELEMENT ENERGY PROJECTION TECHNIQUE[J]. Engineering Mechanics, 2020, 37(9): 1-7, 29. doi: 10.6052/j.issn.1000-4750.2020.07.0446

基于EEP技术的一维有限元结点位移误差计算

doi: 10.6052/j.issn.1000-4750.2020.07.0446
基金项目: 国家自然科学基金项目(51878383,51378293)
详细信息
    作者简介:

    邢沁妍(1981−),女,辽宁人,副教授,博士,主要从事结构工程研究(E-mail: xingqy@tsinghua.edu.cn)

    袁 全(1993−),男,北京人,博士生,主要从事结构工程研究(E-mail: yuanq19@mails.tsinghua.edu.cn)

    通讯作者: 袁 驷(1953−),男,北京人,教授,博士,主要从事结构工程研究(E-mail: yuans@tsinghua.edu.cn)

CALCULATION OF ERRORS OF NODAL DISPLACEMENTS IN ONE-DIMENSIONAL FINITE ELEMENT METHODS USING ELEMENT ENERGY PROJECTION TECHNIQUE

  • 摘要: 利用单元能量投影(Element Energy Projection,简称EEP)法所计算的EEP超收敛解,在不改变有限元网格及其整体刚度矩阵的情况下,导出残差的等效结点荷载向量,只经回代过程即可得到具有更高阶精度的结点位移的误差估计,使结点位移精度得到极大提高。该文以一般的二阶常微分方程边值和初值问题为例,给出算法和相应的数值算例。从中可以看出,本法十分简单而高效:对于$m \geqslant {{1}} $次单元,采用EEP简约格式和凝聚格式修正后的结点位移,分别具有$ O({h^{2m\, + \, 2}})$$ O({h^{3m + od {\kern 1pt} (m, \, \, \, 2)}}) $的超常规的超收敛阶。该文给出了典型算例,并对该法的进一步拓展和应用作了讨论。
  • 图  1  例3的精确解

    Figure  1.  Exact solution of Example 3

    图  2  有限元解误差

    Figure  2.  Errors of FE solution

    图  3  整体修正误差

    Figure  3.  Errors of globally updated solution

    图  4  逐单元修正误差

    Figure  4.  Errors of element-wise updated solution

    表  1  结点收敛阶(简约格式)

    Table  1.   Nodal convergence orders (Simplified form)

    单元数${N_e}$$m = 1$$m = 2$$m = 3$$m = 4$
    $\tilde e_{i,\max }^h$$\rho $$\tilde e_{i,\max }^h$$\rho $$\tilde e_{i,\max }^h$$\rho $$\tilde e_{i,\max }^h$$\rho $
    28.5070$ \times {10^{ - 5}}$4.582.0994$ \times {10^{ - 7}}$6.126.7951$ \times {10^{ - 10}}$8.069.4085$ \times {10^{ - 13}}$10.07
    44.7963$ \times {10^{ - 6}}$4.153.1677$ \times {10^{ - 9}}$6.052.6301$ \times {10^{ - 12}}$8.019.0730$ \times {10^{ - 16}}$10.02
    82.9473$ \times {10^{ - 7}}$4.024.9025$ \times {10^{ - 11}}$6.011.0251$ \times {10^{ - 14}}$8.008.8324$ \times {10^{ - 19}}$10.00
    161.8344$ \times {10^{ - 8}}$4.017.6416$ \times {10^{ - 12}}$6.004.0019$ \times {10^{ - 17}}$8.008.6185$ \times {10^{ - 22}}$10.00
    321.1453$ \times {10^{ - 9}}$4.001.1932$ \times {10^{ - 14}}$6.001.5630$ \times {10^{ - 19}}$8.008.4149$ \times {10^{ - 25}}$10.00
    $\bar \rho \approx $46810
    下载: 导出CSV

    表  2  结点收敛阶(凝聚格式)

    Table  2.   Nodal convergence orders (Condensed form)

    单元数${N_e}$$m = 2$$m = 3$$m = 4$$m = 5$
    $\tilde e_{i,\max }^h$ρ$\tilde e_{i,\max }^h$ρ$\tilde e_{i,\max }^h$ρ$\tilde e_{i,\max }^h$ρ
    29.4026$ \times {10^{ - 8}}$7.151.6605$ \times {10^{ - 11}}$10.133.8765$ \times {10^{ - 11}}$12.549.8103$ \times {10^{ - 20}}$16.24
    41.2064$ \times {10^{ - 9}}$6.281.5218$ \times {10^{ - 14}}$10.098.5672$ \times {10^{ - 14}}$12.141.4350$ \times {10^{ - 24}}$16.06
    81.7893$ \times {10^{ - 11}}$6.071.4623$ \times {10^{ - 16}}$10.022.0392$ \times {10^{ - 16}}$12.042.1667$ \times {10^{ - 29}}$16.01
    162.7591$ \times {10^{ - 13}}$6.011.4223$ \times {10^{ - 19}}$10.014.9470$ \times {10^{ - 19}}$12.013.2973$ \times {10^{ - 34}}$16.00
    324.3152$ \times {10^{ - 15}}$6.001.3874$ \times {10^{ - 22}}$10.001.2058$ \times {10^{ - 22}}$12.005.0279$ \times {10^{ - 39}}$16.00
    $\bar \rho \approx $6101216
    下载: 导出CSV

    表  3  反复修正结点位移

    Table  3.   Repeated nodal value updating

    结点1有限元解修正的结点值结点误差
    $u_{\rm{1}}^h$ 0.2727272727 0.2727272727 0.0116050160
    $\varepsilon _{\rm{1}}^h$ 0.0136363636 0.2863636363 −0.0020313477
    $\delta \varepsilon _{\rm{1}}^h$ −0.0021467926 0.2842168437 0.0001154449
    ${\delta ^2}\varepsilon _{\rm{1}}^h$ 0.0000908271 0.2843076709 0.0000246178
    ${\delta ^{\rm{3}}}\varepsilon _{\rm{1}}^h$ 0.0000300643 0.2843377351 −0.0000054464
    ${\delta ^{\rm{4}}}\varepsilon _{\rm{1}}^h$ −0.0000055624 0.2843321727 0.0000001159
    ··· ··· ··· ···
    精确解 0.2843322887
    下载: 导出CSV

    表  4  4个线性元的结点误差

    Table  4.   Errors of nodal values of four linear elements

    网格${N_e} = 4$,$m = 1$
    结点 $i = 1$ $i = 2$ $i = 3$ $i = 4$
    $\varepsilon _i^h$ −0.1802$ \times {10^{ - 3}}$ −0.2242$ \times {10^{ - 3}}$ 0.2363$ \times {10^{ - 4}}$ 0.6194$ \times {10^{ - 3}}$
    $\varepsilon _i^h + \delta \varepsilon _i^h$ −0.1776$ \times {10^{ - 3}}$ −0.2198$ \times {10^{ - 3}}$ 0.2017$ \times {10^{ - 4}}$ 0.6153$ \times {10^{ - 3}}$
    $e_i^h$ −0.1775$ \times {10^{ - 3}}$ −0.2197$ \times {10^{ - 3}}$ 0.2022$ \times {10^{ - 4}}$ 0.6146$ \times {10^{ - 3}}$
    下载: 导出CSV

    表  5  结点误差和收敛阶($m = 1$,简约格式)

    Table  5.   Nodal errors and convergence orders ($m = 1$, simplified form)

    单元步长$h$有限元解$u_i^h$整体修正解$u_i^h + \varepsilon _i^h$逐单元修正解${(u_i^h + \varepsilon _i^h)^*}$
    $\mathop {\max }\limits_i \left| {u_i^{} - u_i^h} \right|$$\rho $$\mathop {\max }\limits_i \left| {u_i^{} - u_i^h - \varepsilon _i^h} \right|$$\rho $$\mathop {\max }\limits_i \left| {u_i^{} - {{(u_i^h + \varepsilon _i^h)}^*}} \right|$$\rho $
    0.496.878$ \times {10^{ - 3}}$23.437$ \times {10^{ - 3}}$1.449$ \times {10^{ - 3}}$
    0.224.400$ \times {10^{ - 3}}$1.991.484$ \times {10^{ - 3}}$3.980.095$ \times {10^{ - 3}}$3.92
    0.16.102$ \times {10^{ - 3}}$2.000.093$ \times {10^{ - 3}}$4.000.006$ \times {10^{ - 3}}$3.98
    $\bar \rho \approx $244
    下载: 导出CSV
  • [1] Strang G, Fix G J. An analysis of the finite element method [M]. New Jersey: Prentice-Hall, 1973.
    [2] Douglas J, Dupont T. Galerkin approximations for the two point boundary problems using continuous piecewise polynomial spaces [J]. Numerische Mathematik, 1974, 22(2): 99 − 109. doi:  10.1007/BF01436724
    [3] Hulme B L. One-step piecewise polynomial galerkin methods for initial value problems [J]. Mathematics of Computation, 1972, 26(118): 415 − 426. doi:  10.1090/S0025-5718-1972-0321301-2
    [4] Zhu J Z, Zienkiewicz O C. Superconvergence recovery technique and a posteriori error estimator [J]. International Journal for Numerical Methods in Engineering, 1990(30): 1321 − 1339.
    [5] Zienkiewicz O C, Zhu J Z. The superconvergence patch recovery and a posteriori error estimator, part I: the superconvergence patch recovery [J]. International Journal for Numerical Methods in Engineering, 1992(33): 1331 − 1364.
    [6] 袁驷. 从矩阵位移法看有限元应力精度的损失于恢复[J]. 力学与实践, 1998, 20(4): 1 − 6.

    Yuan Si. The loss and recovery of stress accuracy in FEM as seen from matrix displacement method [J]. Mechanics in Engineering, 1998, 20(4): 1 − 6. (in Chinese)
    [7] 袁驷, 王枚. 一维有限元后处理超收敛解答计算的EEP法[J]. 工程力学, 2004, 21(2): 1 − 9. doi:  10.3969/j.issn.1000-4750.2004.02.001

    Yuan Si, Wang Mei. An element-energy-projection method for post-computation of super-convergent solutions in one-dimensional FEM [J]. Engineering Mechanics, 2004, 21(2): 1 − 9. (in Chinese) doi:  10.3969/j.issn.1000-4750.2004.02.001
    [8] 袁驷, 林永静. 二阶非自伴两点边值问题Galerkin有限元后处理超收敛解答计算的EEP法[J]. 计算力学学报, 2007, 24(2): 142 − 147. doi:  10.3969/j.issn.1007-4708.2007.02.004

    Yuan Si, Lin Yongjing. An EEP method for post-computation of super-convergent solutions in one-dimensional Galerkin FEM for second order non-self-adjont BVP [J]. Journal of Computational Mechanics, 2007, 24(2): 142 − 147. (in Chinese) doi:  10.3969/j.issn.1007-4708.2007.02.004
    [9] 袁驷, 王旭, 邢沁妍, 叶康生. 具有最佳超收敛阶的EEP法计算格式:Ⅰ算法公式[J]. 工程力学, 2007, 24(10): 1 − 5. doi:  10.3969/j.issn.1000-4750.2007.10.001

    Yuan Si, Wang Xu, Xing Qinyan, Ye Kangsheng. A scheme with optimal order of super-convergence based on the element energy projection method-I Formulation [J]. Engineering Mechanics, 2007, 24(10): 1 − 5. (in Chinese) doi:  10.3969/j.issn.1000-4750.2007.10.001
    [10] Yuan S, Wu Y, Xu J J, Yuan Z, Xing Q Y, Ye K S. A super-convergence strategy for two-dimensional FEM based on element energy projection technique [J]. Journal of Nanoelectronics and Optoelectronics, 2017, 12(11): 1284 − 1294. doi:  10.1166/jno.2017.2272
    [11] Yuan S, Wu Y, Xing Q Y. Recursive super-convergence computation for multi-dimensional problems via one-dimensional element energy projection technique [J]. Applied Mathematics and Mechanics, 2018, 39(7): 1031 − 1044. doi:  10.1007/s10483-018-2345-7
    [12] Yuan S, Ye K S, Wang Y L, Kennedy D, Williams F W. Adaptive finite element method for eigensolutions of regular second and fourth order Sturm-Liouville problems via the element energy projection technique [J]. Engineering Computations, 2017, 34(8): 2862 − 2876. doi:  10.1108/EC-03-2017-0090
    [13] Dong Y Y, Yuan S, Xing Q Y. Adaptive finite element analysis with local mesh refinement based on a posteriori error estimate of element energy projection technique [J]. Engineering Computations, 2019, 36(6): 2010 − 2033.
    [14] Zhao Q H, Zhou S Z, Zhu Q D. Mathematical analysis of EEP method for one-dimensional finite element postprocessing [J]. Applied Mathematics and Mechanics, 2007, 28(4): 441 − 445. doi:  10.1007/s10483-007-0403-y
    [15] 袁驷, 赵庆华. 具有最佳超收敛阶的EEP法计算格式: Ⅲ数学证明[J]. 工程力学, 2007, 24(12): 1 − 6. doi:  10.3969/j.issn.1000-4750.2007.12.001

    Yuan Si, Zhao Qinghua. A scheme with optimal order of super-convergence based on the element energy projection method-III Mathematical analysis [J]. Engineering Mechanics, 2007, 24(12): 1 − 6. (in Chinese) doi:  10.3969/j.issn.1000-4750.2007.12.001
  • [1] 孙浩涵, 袁驷.  基于EEP超收敛解的自适应有限元法特性分析 . 工程力学, 2019, 36(2): 17-25. doi: 10.6052/j.issn.1000-4750.2018.10.0439
    [2] 叶康生, 姚葛亮.  平面曲梁有限元静力分析的p型超收敛算法 . 工程力学, 2017, 34(11): 26-33,58. doi: 10.6052/j.issn.1000-4750.2016.07.0572
    [3] 钟博, 叶康生, 袁驷.  一维C1有限元后处理超收敛计算的新型算法 . 工程力学, 2017, 34(9): 27-33, 53. doi: 10.6052/j.issn.1000-4750.2016.05.0363
    [4] 叶康生, 曾强.  结构自由振动问题有限元新型超收敛算法研究 . 工程力学, 2017, 34(1): 45-50,68. doi: 10.6052/j.issn.1000-4750.2015.05.0421
    [5] 徐俊杰, 袁驷.  有限元线法EEP超收敛计算简约格式的再简约 . 工程力学, 2016, 33(增刊): 1-5. doi: 10.6052/j.issn.1000-4750.2015.05.S025
    [6] 钟博, 叶康生, 袁驷.  基于p型超收敛计算的一维有限元自适应分析 . 工程力学, 2016, 33(增刊): 23-28. doi: 10.6052/j.issn.1000-4750.2015.05.S038
    [7] 袁驷, 吴越, 徐俊杰, 邢沁妍.  基于EEP法的三维有限元超收敛计算初探 . 工程力学, 2016, 33(9): 15-20. doi: 10.6052/j.issn.1000-4750.2016.05.ST07
    [8] 袁驷, 刘泽洲, 邢沁妍.  一维变分不等式问题的自适应有限元分析新探 . 工程力学, 2015, 32(7): 11-16. doi: 10.6052/j.issn.1000-4750.2014.07.ST05
    [9] 袁驷, 邢沁妍, 叶康生.  一维C1有限元EEP超收敛位移计算简约格式的误差估计 . 工程力学, 2015, 32(9): 16-19. doi: 10.6052/j.issn.1000-4750.2015.03.0695
    [10] 袁驷, 邢沁妍.  一维Ritz有限元超收敛计算的EEP法简约格式的误差估计 . 工程力学, 2014, 31(12): 1-3,16. doi: 10.6052/j.issn.1000-4750.2014.11.0973
    [11] 唐义军, 罗建辉.  一维有限单元法位移模式研究 . 工程力学, 2013, 30(3): 46-52,57. doi: 10.6052/j.issn.1000-4750.2011.09.0629
    [12] 袁 驷, 肖 嘉, 叶康生.  线法二阶常微分方程组有限元分析的EEP超收敛计算 . 工程力学, 2009, 26(11): 1-009,.
    [13] 赵庆华, 周叔子.  关于单元能量投影法的两点注记 . 工程力学, 2008, 25(2): 0-094,.
    [14] 袁 驷, 邢沁妍, 叶康生.  具有最佳超收敛阶的Galerkin有限元EEP法计算格式 . 工程力学, 2008, 25(11): 1-007.
    [15] 袁 驷, 王 枚, 王 旭.  二维有限元线法超收敛解答计算的EEP法 . 工程力学, 2007, 24(1): 0-010.
    [16] 袁 驷, 邢沁妍, 王 旭, 叶康生.  具有最佳超收敛阶的EEP法计算格式:II 数值算例 . 工程力学, 2007, 24(11): 0-006.
    [17] 袁 驷, 赵庆华.  具有最佳超收敛阶的EEP法计算格式: III 数学证明 . 工程力学, 2007, 24(12): 0-005,.
    [18] 袁 驷, 王 旭, 邢沁妍, 叶康生.  具有最佳超收敛阶的EEP法计算格式: I 算法公式 . 工程力学, 2007, 24(10): 0-005.
    [19] 袁驷, 王枚, 和雪峰.  一维C1有限元超收敛解答计算的EEP法 . 工程力学, 2006, 23(2): 1-9.
    [20] 袁驷, 王枚.  一维有限元后处理超收敛解答计算的EEP法 . 工程力学, 2004, 21(2): 1-9.
  • 加载中
图(4) / 表 (5)
计量
  • 文章访问数:  56
  • HTML全文浏览量:  20
  • PDF下载量:  31
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-07-14
  • 修回日期:  2020-07-21
  • 网络出版日期:  2020-09-07
  • 刊出日期:  2020-09-25

基于EEP技术的一维有限元结点位移误差计算

doi: 10.6052/j.issn.1000-4750.2020.07.0446
    基金项目:  国家自然科学基金项目(51878383,51378293)
    作者简介:

    邢沁妍(1981−),女,辽宁人,副教授,博士,主要从事结构工程研究(E-mail: xingqy@tsinghua.edu.cn)

    袁 全(1993−),男,北京人,博士生,主要从事结构工程研究(E-mail: yuanq19@mails.tsinghua.edu.cn)

    通讯作者: 袁 驷(1953−),男,北京人,教授,博士,主要从事结构工程研究(E-mail: yuans@tsinghua.edu.cn)

摘要: 利用单元能量投影(Element Energy Projection,简称EEP)法所计算的EEP超收敛解,在不改变有限元网格及其整体刚度矩阵的情况下,导出残差的等效结点荷载向量,只经回代过程即可得到具有更高阶精度的结点位移的误差估计,使结点位移精度得到极大提高。该文以一般的二阶常微分方程边值和初值问题为例,给出算法和相应的数值算例。从中可以看出,本法十分简单而高效:对于$m \geqslant {{1}} $次单元,采用EEP简约格式和凝聚格式修正后的结点位移,分别具有$ O({h^{2m\, + \, 2}})$$ O({h^{3m + od {\kern 1pt} (m, \, \, \, 2)}}) $的超常规的超收敛阶。该文给出了典型算例,并对该法的进一步拓展和应用作了讨论。

English Abstract

袁驷, 邢沁妍, 袁全. 基于EEP技术的一维有限元结点位移误差计算[J]. 工程力学, 2020, 37(9): 1-7, 29. doi: 10.6052/j.issn.1000-4750.2020.07.0446
引用本文: 袁驷, 邢沁妍, 袁全. 基于EEP技术的一维有限元结点位移误差计算[J]. 工程力学, 2020, 37(9): 1-7, 29. doi: 10.6052/j.issn.1000-4750.2020.07.0446
Si YUAN, Qin-yan XING, Quan YUAN. CALCULATION OF ERRORS OF NODAL DISPLACEMENTS IN ONE-DIMENSIONAL FINITE ELEMENT METHODS USING ELEMENT ENERGY PROJECTION TECHNIQUE[J]. Engineering Mechanics, 2020, 37(9): 1-7, 29. doi: 10.6052/j.issn.1000-4750.2020.07.0446
Citation: Si YUAN, Qin-yan XING, Quan YUAN. CALCULATION OF ERRORS OF NODAL DISPLACEMENTS IN ONE-DIMENSIONAL FINITE ELEMENT METHODS USING ELEMENT ENERGY PROJECTION TECHNIQUE[J]. Engineering Mechanics, 2020, 37(9): 1-7, 29. doi: 10.6052/j.issn.1000-4750.2020.07.0446
  • 有限元后处理超收敛计算是近年来研究热点之一[1-5]。根据一维有限元的超收敛理论,无论是边值问题的Ritz有限元还是初值问题的Galerkin有限元的位移解,$m\;( \geqslant {\rm{1)}}$次多项式单元可得到具有${h^{2m}}$阶的超收敛端结点位移(真解足够光滑),然而单元内部的有限元解仅为${h^{m + 1}}$[1-3]。袁驷等[6-9]提出了有限元后处理中超收敛计算的单元能量投影(Element Energy Projection,简称EEP)法,能够计算一维有限元的逐点超收敛解。该技术虽然是基于一维有限元的,但通过一种逐维(dimension-by-dimension,简称D-by-D)恢复策略[10-11],已经拓展到二维和三维的有限元分析,并应用于多种问题的自适应求解[12-13]

    一维问题的EEP公式有两种格式:简约格式和凝聚格式[7-9]。二者均被证明可以提供单元内任意点位移及其导数的超收敛解。其中,简约格式的超收敛解精度[14]:位移为$O({h^{m + 2}})$ $(m > 1)$,导数为$O({h^{m + 1}})$;凝聚格式的精度[15]:位移及导数均为$O({h^{2m}})$。EEP解并不改变单元端结点的位移值,其精度已经很高了,即$O({h^{2m}})$。然而,更高精度的结点位移对于结点误差估计及其精度提升均具有重要意义。例如,对于初值问题,结点位移精度对于稳定可靠的时程积分是至关重要的,因此,需要有能够对结点位移误差进行有效估计和控制的技术。

    本文充分利用EEP解超收敛特性,提出一种简便有效的结点位移误差计算方法。本法在传统有限元计算后,在不改变现有网格和整体刚度矩阵的情况下,用EEP解的残余荷载生成新的荷载向量,只经常规的回代过程,即可得到高精度的结点位移误差。再将所估计的误差叠加在传统有限元解,会得到修正的结点位移,其精度具有更高阶的超收敛性:对简约格式有$O({h^{2m + 2}})$精度,对凝聚格式有$O({h^{3m + od \; (m,\,2)}})$精度。例如,对于线性元,修正后的结点位移具有$O({h^4})$精度。此外,由于修正后的结点位移再次具有更高阶的超收敛性,因此可以用其再次改进EEP解,进而计算另一轮超收敛结点位移。如此反复迭代,最终可获得精确的结点位移,其中每个迭代步的计算量只是荷载向量的生成和回代。

    本文以一般的二阶常微分方程为模型问题,给出本法的基本思路和算法,并给出典型算例以展示本法优越的超收敛性以及多样化的拓展和应用。

    • 本文考虑一般的二阶常微分方程边值问题(BVP)和初值问题(IVP):

      $$ \left\{ \begin{array}{l} Lu \equiv - (p(x)u')' + r(x)u' + q(x)u = f(x),\\ \quad \quad \quad p(x) > 0,\;\quad 0 < x < 1\\ {\rm{BVP:}}\;\; u(0) = 0,\quad u'(1) = 0\\ {\rm{IVP:}}\;\;\;u(0) = {u_0},\;\; u'({\rm{0}}) = {v_0} \end{array} \right.\;\;\;\;\;\;\;\begin{array}{*{20}{c}} {\left( {1{\rm{a}}} \right)}\\ {\left( {1{\rm{b}}} \right)} \end{array} $$

      式中:$(\;)' = {{{\rm{d(}}\,{\rm{)}}} / {{\rm{d}}x}}$$L$为式(1)中定义的线性微分算子。不失一般性,以下将$u(x)$称为位移。对于非自伴情况,其伴随算子为:

      $$ {L^*}v \equiv - (p(x)v')' - (r(x)v)' + q(x)v $$ (2)

      定义双线性型和线性型如下:

      $$ \begin{split} & a(u,v) = \int_{\,0}^{\,1} {\;(pu'v' + ru'v + quv)\,{\rm{d}}x}, \\& (f,v) = \int_{\,0}^{\,1} { fv\,{\rm{d}}x} - p(0){v_0}v(0) \end{split} $$ (3)

      则用常规的Galerkin位移型有限元求解上述问题归结为求解如下的弱形式问题:

      $$ a({u^h},{v^h}) = (f,{v^h}),\;\;\;\;\forall {v^h} \in {S^h} $$ (4)

      式中:${u^h} \in {S^h}$为有限元解;${S^h}$m次分段多项式所组成的有限元试探(检验)空间。由${v^h}$的任意性,可导出通常意义下的刚度方程:

      $$ {{K}}\, {{D = }}\,{{P}} $$ (5)

      式中:${{K}} $为整体刚度矩阵;${{D}}$为结点位移向量;${{P}}$为等效结点荷载向量。为方便起见,本文中的“结点”在未特别提及时均指“单元两端结点”,而相应的结点位移表示为$u_i^h$

      根据一维有限元误差估计理论,当问题式(1)的解足够光滑且网格足够密(单元长度h足够小)时,由m次元(试探函数和检验函数同次)得到的有限元解在单元端结点上获得$O({h^{2m}})$[2-3]的超收敛性,而在单元内部则为常规的$O({h^{m + 1}})$[1]的收敛性,即:

      $$ \left| {{u_i} - u_i^h} \right| \leqslant C{h^{2m}}{\left\| u \right\|_{m + 1}}\;\;\;\; $$ (6)
      $$ {\left\| {u - {u^h}} \right\|_0} \leqslant C{h^{m + 1}}{\left\| u \right\|_{m + 1}} $$ (7)

      式中,$|| \cdot |{|_0}$$|| \cdot |{|_{m + 1}}$分别为整个解域上${H^0}$${H^{m + 1}}$范数。式(7)的超收敛阶已经非常高,通常认为是最佳超收敛阶。本文提出用同次单元、同一网格进行另一轮有限元分析而计算出结点位移误差${u_i} - u_i^h$近似解的方法,从而可获得精度更高的结点位移修正解。这一方法关键技术是EEP超收敛解的应用,将在下一小节中讨论。

    • 袁驷等[6]基于结构力学中的矩阵位移法和有限元中的投影定理[1],提出了一种一维有限元后处理超收敛计算的新型方法,简称为EEP法[7-9]。本节简要介绍该方法的基本思想和相关公式。

      记模型问题式(1)的常规有限元解的误差为${e^h} = u - {u^h}$。由有限元数学理论中熟知的投影定理有$a({e^h},{v^h}) = {\rm{0}}$,$\forall {v^h} \in {S^h}$。值得注意的是,投影定理是在整个解域$[{\rm{0}},{\rm{1}}]$上定义的,而不是在单个单元上定义的。然而,EEP方法的主要思想是假设投影定理在单个单元$e = [{\bar x_1},{\bar x_2}]$ (${\bar x_1}$${\bar x_2}$为两端点坐标)上几乎成立,进而数学分析证明了单元投影的误差为$O({h^{2m}})$[10],即:

      $$ \begin{split} & {a^e}({e^h},{v^h}) = \int_{\;{{\bar x}_1}}^{\;{{\bar x}_2}} {(p{e^{h\prime}} {v^{h\prime}} + r{e^{h\prime}} {v^h} + q{e^h}{v^h}{\rm{)d}}x} \leqslant C{h^{2m}}, \\ & \forall {v^h} \in S_e^h = \{ {v^h}| {{v^h} = \sum\nolimits_{i = 1}^{m + 1} {{N_i}} v_i^h,v_i^h \in R} \} \\[-16pt] \end{split} $$ (8)

      式中,${N_i}$($i = 1,\;2,\cdots,m + 1$)为m次拉格朗日形函数或升阶谱形函数。由整体投影定理可知,对于内部形函数有${a^e}({e^h},{N_j}) = 0$ ($j = 2,3,\cdots,m$),则式(8)也可仅用线性检验函数来表示如下:

      $$ \begin{split} & {a^e}({e^h},{v^h}) \leqslant C{h^{2m}} ,\\ & \forall {v^h} \in S_e^h = \{ {v^h}| {{v^h} = \sum\nolimits_{i = 1}^2 {{{\bar N}_i}} v_i^h,v_i^h \in R} \} \end{split} $$ (9)

      式中,${\bar N_i}$($i = 1,\;2$)为线性形函数。单元能量投影(EEP)的名称就是从${a^e}({e^h},{v^h}) \approx 0$的假设而来,而这一假设在各类一维有限元中已被证明其合理性。

      EEP具有两种位移超收敛计算格式[7-9]

      1) 简约格式:

      $$\tag{10a}\begin{split} u_a^* = & u_a^h + \frac{h}{{{p_a}}}\left( {{{\bar N}_{1a}}\int_{\,{{\bar x}_1}}^{\,{{\bar x}_a}} {(f - L{u^h}){{\bar N}_2}{\rm{d}}x} +} \right. \\ & \left. { {{\bar N}_{2a}}\int_{\,{{\bar x}_a}}^{\,{{\bar x}_2}} {(f - L{u^h}){{\bar N}_1}{\rm{d}}x} } \right) \end{split} $$ (10a)

      2) 凝聚格式:

      $$\tag{10b} \begin{split} u_a^* =& u_a^h + \frac{1}{{{p_a}{{\tilde v}_a}}}\left( {\tilde N_{1a}^*\int_{\;{{\bar x}_1}}^{\;{{\bar x}_a}} {(f - L{u^h})\,\tilde N_2^*{\rm{d}}x} +} \right. \\ & \left. { \tilde N_{2a}^*\int_{\;{{\bar x}_a}}^{\;{{\bar x}_2}} {(f - L{u^h})\,\tilde N_1^*{\rm{d}}x} } \right) ,\\ {\tilde v_a} =& \tilde N_{1a}^*{(\tilde N_2^*)'_a} - \tilde N_{2a}^*{(\tilde N_1^*)'_a} \end{split} $$ (10b)

      式中:$h = {\bar x_2} - {\bar x_1}$为单元长度;${(\;)^{*}}$${(\;)^h}$${(\;)_a}$分别为超收敛解、有限元解和在${\bar x_a} \in [{\bar x_1},{\bar x_2}]$点的值;$\tilde N_i^*$ ($i = 1,\;2$)为两个m次凝聚形函数。凝聚形函数是通过“静力凝聚”手段将内部结点凝聚后获得,也等效于用单个m次单元对精确形函数$N_{Ei}^*$ ($i = 1,\;2$) (精确形函数满足:${L^*}N_{Ei}^* = 0$$N_{Ei}^*({\bar x_j}) = {\delta _{ij}}$ $(j = 1,2)$)求得的有限元近似解。可以看出,如果用${\bar N_i}$代替式(10b)的${\tilde N_i}$即为式(10a)。更多凝聚形函数内容请参考文献[9]。EEP导数公式也可参见文献[7-9]。

      对于上述EEP解,已经得到数学理论证明,即对于m次单元,用式(10a)和式(10b)计算出的超收敛位移有以下误差估计[14-15]

      $$ \begin{split} & {\text{简约格式}}:\;\left\{\!\! {\begin{array}{*{20}{l}} {\left| {u - u_{}^*} \right| \leqslant C{h^{\rm{2}}}||u|{|_{\rm{2}}}},\quad\quad\;\;\;\; m=1 \\ {\left| {u - u_{}^*} \right| \leqslant {C_1}{h^{m + {\rm{2}}}}||u|{|_{m + 1}}},\;\; m > 1 \end{array}} \right.\\& {\text{凝聚格式}}:\;\left| {u - u_{}^*} \right| \leqslant C{h^{2m}}||u|{|_{m + 1}},\;\;\quad \quad m \geqslant 1 \end{split}\!\!\!\!\!\! $$ (11)

      式中,$C$${C_1}$表示通常的常数。可见,对于简约格式,EEP解${u^*}$在精度上至少比传统有限元解${u^h}$高一阶;而对于凝聚格式,EEP解${u^*}$与结点位移$u_i^h$同阶。

      此外,上述EEP方案还有一些优秀的性质,其中与本文有关的有:

      a) 当${\bar x_a}$位于端点时(即${\bar x_a} = {\bar x_i}$, $i = 1,\;2$),$u_{}^*$与有限元解$u_{}^h$相同,即$u_i^* = u_i^h$。这意味着有限元解${u^h}$可看作对EEP解${u^*}$在单元端点的插值。换言之,EEP解仅修复了单元内部的精度,因此有限元结点位移误差$e_i^h = {u_i} - u_i^h$等同于EEP解的误差$e_i^* = {u_i} - u_i^*$,即$e_i^* = e_i^h$

      b) 对于公共结点,结点两侧单元计算出的EEP导数$u{_i^{*\prime} }$是连续的。

      c) 特别地,在自由边界点,EEP导数$u{_i^{*\prime} }$满足精确值,亦即EEP解满足自然边界条件。

      d) EEP超收敛计算是一种后处理方法,因此,不依赖于有限元求解过程。

    • 理想情况下,有限元误差是计算${e^h} = u - {u^h}$。一个自然的想法,是将${e^h}$回代到原始微分方程中得到新的边值(初值)问题,即求解:

      $$ \left\{ \begin{aligned} & L{e^h} = {r^h},\quad 0 < x < 1\\& {\rm{BVP:}}\;{e^h}(0) = 0,\quad {e^{h\prime}} (1) = - {u^{h\prime}} (1)\\& {\rm{IVP:}}\;\;\;{e^h}(0) = 0,\quad {e^{h\prime}} ({\rm{0}}) = {v_0} - {u^{h\prime}} ({\rm{0}}) \end{aligned} \right. $$ (12)

      式中:残余荷载${r^h} = f - L{u^h}$。但是,如果不改变原始网格和单元次数,而去求解${e^h}$是徒劳的,只能得到零解。另一方面,若网格加密后再求解${e^h}$,则计算成本将超过求解原问题本身。因此,直接求解${e^h}$不可取。

      注意到有限元解${u^h}$与EEP解${u^*}$在单元端结点处相同,即$u_i^* = u_i^h$,则如果误差${e^*} = u - {u^*}$可以计算,在端点处$e_i^* = e_i^h$。因此,用相同网格、同次单元计算${e^*}$的有限元近似解${\varepsilon ^h}$,其结点值$\varepsilon \; _i^h$就是有限元结点位移误差的近似解,而且由于其解为结点值,仍具有超收敛性。此为本文的基本思想,简单有效,以下给出具体做法。

    • ${e^*} = u - {u^*}$代入问题式(1),可得到${e^*}$的定解方程和边界条件:

      $$ \left\{ \begin{array}{l} L{e^*} = {r^*},\quad 0 < x < 1\\ {\rm{BVP:}}\;\;{e^*}(0) = 0,\quad {e^{*\prime}} (1) = 0\\ {\rm{IVP:}}\;\;\;\;{e^*}(0) = 0,\quad {e^{*\prime}} ({\rm{0}}) = 0 \end{array} \right. $$ (13)

      式中,残差${r^*} = f - L{u^*}$,而各点的${u^*}$可用式(10)算出。相比问题式(1),除了荷载和边界条件不同,其他定解条件均相同。这意味着,在原有网格上求解$e{\; ^*}$的近似解${\varepsilon ^h}$,整体刚度矩阵${{K}}\,$及其逆矩阵${{{K}}^{ - 1}}$都不变,只是等效结点荷载向量${{P}}$需重新生成为${{{P}}^{*}}$,再经一个回代过程即可求得结点位移误差向量$\delta {{D}}$,进而可提取$\varepsilon _i^h$

    • 上述讨论可归纳为两阶段的算法:

      阶段1:常规有限元解

      1)在当前网格建立有限元刚度方程${{K}}\,{{D}} = {{P}}$

      2)解${{D}} = {{{K}}^{ - 1}}\,{{P}}$,得到有限元解${u^h}$

      阶段2:结点位移修正

      3)计算各单元超收敛解${u^*}$,用残差${r^*} = f - L{u^*}$替代$f$生成结点荷载向量${{{P}}^*}$

      4)解$\delta {{D}} = {{{K}}^{ - 1}}\,{{{P}}^*}$(仅回代即可),得有限元解${{D}}$的误差近似解$\delta {{D}}$

      5)更新结点位移${{D}} + \delta {{D}}$$u_i^h + \varepsilon _i^h$

      注意,算法的最后一步意味着本法也是一个结点修正法。在实际计算中,步骤2)一般采用矩阵的三角分解,并不直接计算逆矩阵${{{K}}^{ - 1}}$,故步骤4)无需再作矩阵分解,只需回代过程即可。

    • 第2.3节算法是一个基本算法,应用中可以有多种变化和拓展,简述几点如下:

      1)结点位移全量。记${{{D}}^{(1)}} = {{D}} + \delta {{D}}$,算法步骤4)可改为${{{D}}^{(1)}} = {{{K}}^{ - 1}}\,({{P}} + {{{P}}^*})$,亦即本法不仅可以计算结点位移增量,也能直接计算修正后的结点位移全量。

      2)修正的EEP解。在求解了$\delta {{D}}$得到${\varepsilon ^h}$后,还可以用EEP法再计算全新的超收敛解${\varepsilon ^*}$,此即为${u^*}$的逐点误差估计,从而得到更精确的有限元逐点误差估计${e^h} \approx {u^*} + {\varepsilon ^*} - {u^h}$

      3)结点的再修正。继得到上述${\varepsilon ^*}$后,还可以用本文第2.3节的算法再计算${\varepsilon ^*}$的误差$\delta \; {e^*} = {e^*} - {\varepsilon ^*}$的近似解$\delta {\varepsilon ^h}$(同时也是修正的位移全量的误差近似解)。$\delta \; {e^*}$的定解方程和条件为:

      $$\left\{ \begin{array}{l} L\,\delta {e^*} = \delta {r^*},\quad 0 < x < 1\\ {\rm{BVP:}}\;\;\delta {e^*}(0) = 0,\quad \delta {e^{*\prime}} (1) = 0\\ {\rm{IVP:}}\;\;\;\;\delta {e^*}(0) = 0,\quad \delta {e^{*\prime}} ({\rm{0}}) = 0 \end{array} \right.$$ (14)

      式中,残差$\delta {r^*} = {r^*} - L{\varepsilon ^*}$。用“荷载”$\delta {r^*}$生成其等效荷载向量$\delta {{{P}}^*}$,无需重新“求逆”,可求得更高阶的结点位移误差的近似解${\delta ^2}{{D}} = {{{K}}^{ - 1}}\,\delta {{P}}{\; ^*}$,也记作$\delta \varepsilon _i^h$

      4)解的反复修正。继续上述迭代和叠加,若迭代过程收敛,相当于可以通过下面矩阵方程得到精确结点位移向量${{{D}}_{\rm{E}}}$

      $$ {{K}}{{{D}}_{\rm{E}}} = {{{P}}_{\rm{I}}} $$

      其中:

      $$ {{{P}}_{\rm{I}}} = {{P}} + {{{P}}^*} + \delta {{{P}}^*} + {\delta ^2}{{{P}}^*} + \cdots $$ (16)

      如果荷载项依次包含上式右端项,则最终会得到精确的结点位移及其误差:

      $$\tag{17a}{{{D}}_{\rm{E}}} = {{D}} + \delta {{D}} + {\delta ^2}{{D}} + \cdots $$
      $$\tag{17b} e_i^h = \varepsilon _i^h + \delta \varepsilon _i^h + {\delta ^2}\varepsilon _i^h + \cdots $$

      亦即,在不改变原有的单元网格和次数,从而不改变整体刚度矩阵的情况下,可以通过逐步迭代构造一个适当的荷载向量${{{P}}_{\rm{I}}}$,其对应的解是精确的结点位移${{{D}}_{\rm{E}}}$

    • 本节给出若干典型算例,用以验证本法优良的内在特性和多样化的应用场景。本文采用符号计算软件MAPLE计算,其优点在于用户可以定义数值的字长。

      例1. 收敛阶

      记修正的结点解为$\tilde u_i^h = u_i^* + \varepsilon _i^h = u_i^h + \varepsilon _i^h$。对于足够光滑的解答,尽管尚未有严格的数学证明,大量数值算例一致验证了以下的收敛阶:

      $$\tag{18a} {\text{简约格式}}:\;\tilde e_i^h \equiv \left| {{u_i} - \tilde u_i^h} \right| \leqslant C{h^{2m + 2}},\qquad\qquad $$
      $$\tag{18b} {\text{凝聚格式}}:\;\tilde e_i^h \equiv \left| {{u_i} - \tilde u_i^h} \right| \leqslant C{h^{{\rm{3}}m + od \,(m,\,2)}},\;m \geqslant 1 $$

      下面给出数值算例来验证以上的收敛阶。

      计算中,采用${N_e}$ (=1, 2, 4, ···)个单元的均分网格,结点依次为$0 = {x_0} < {x_1} < {x_2} < \cdots < {x_{{N_e}}} = 1$,则单元长度为$h = {1 / {{N_e}}}$。结点位移误差表示为:

      $$ \tilde e_{i,\max }^h \equiv \mathop {\max }\limits_{1 \leqslant i \leqslant {N_e}} \left| {{u_i} - \tilde u_i^h} \right| = \mathop {\max }\limits_{1 \leqslant i \leqslant {N_e}} \left| {e_i^* - \varepsilon _i^h} \right| $$ (19)

      $\rho $${h^\rho }$项中的收敛阶,其估算方法是通过单元长度h均匀减半(二分加密)后两个网格上的结点误差比较推算出。相应地, $\bar \rho $表示各个网格估计的$\rho $的收敛值。

      为简便,取如下的典型边值问题为例:

      $$ \left\{ \begin{array}{l} - \,u'' + u' + \; \; u = 1, \quad \quad 0 < x < 1 \\ u(0) = 0,\quad \quad\quad\quad\quad\quad u'(1) = 0 \\ \end{array} \right. $$ (20)

      其精确解为$u = {{({\alpha _1}{e^{{\alpha _1} + {\alpha _2}x}} - {\alpha _2}{e^{{\alpha _1}x + {\alpha _2}}})} / \beta } + 1$,其中:${\alpha _1},{\,_2} = {{({\rm{1}} \pm \sqrt 5 )} / 2}$$\beta = {\alpha _2}{e^{{\alpha _2}}} - {\alpha _1}{e^{{\alpha _1}}}$。需要说明的是,变系数问题也具有相同的收敛阶。$m = 1\sim 4$次单元的计算结果见表1(简约格式)和表2(凝聚格式),其结果验证了式(18)的有效性。值得注意的是,EEP凝聚格式在$m$为奇数时,会产生附加的收敛阶。

      表 1  结点收敛阶(简约格式)

      Table 1.  Nodal convergence orders (Simplified form)

      单元数${N_e}$$m = 1$$m = 2$$m = 3$$m = 4$
      $\tilde e_{i,\max }^h$$\rho $$\tilde e_{i,\max }^h$$\rho $$\tilde e_{i,\max }^h$$\rho $$\tilde e_{i,\max }^h$$\rho $
      28.5070$ \times {10^{ - 5}}$4.582.0994$ \times {10^{ - 7}}$6.126.7951$ \times {10^{ - 10}}$8.069.4085$ \times {10^{ - 13}}$10.07
      44.7963$ \times {10^{ - 6}}$4.153.1677$ \times {10^{ - 9}}$6.052.6301$ \times {10^{ - 12}}$8.019.0730$ \times {10^{ - 16}}$10.02
      82.9473$ \times {10^{ - 7}}$4.024.9025$ \times {10^{ - 11}}$6.011.0251$ \times {10^{ - 14}}$8.008.8324$ \times {10^{ - 19}}$10.00
      161.8344$ \times {10^{ - 8}}$4.017.6416$ \times {10^{ - 12}}$6.004.0019$ \times {10^{ - 17}}$8.008.6185$ \times {10^{ - 22}}$10.00
      321.1453$ \times {10^{ - 9}}$4.001.1932$ \times {10^{ - 14}}$6.001.5630$ \times {10^{ - 19}}$8.008.4149$ \times {10^{ - 25}}$10.00
      $\bar \rho \approx $46810

      表 2  结点收敛阶(凝聚格式)

      Table 2.  Nodal convergence orders (Condensed form)

      单元数${N_e}$$m = 2$$m = 3$$m = 4$$m = 5$
      $\tilde e_{i,\max }^h$ρ$\tilde e_{i,\max }^h$ρ$\tilde e_{i,\max }^h$ρ$\tilde e_{i,\max }^h$ρ
      29.4026$ \times {10^{ - 8}}$7.151.6605$ \times {10^{ - 11}}$10.133.8765$ \times {10^{ - 11}}$12.549.8103$ \times {10^{ - 20}}$16.24
      41.2064$ \times {10^{ - 9}}$6.281.5218$ \times {10^{ - 14}}$10.098.5672$ \times {10^{ - 14}}$12.141.4350$ \times {10^{ - 24}}$16.06
      81.7893$ \times {10^{ - 11}}$6.071.4623$ \times {10^{ - 16}}$10.022.0392$ \times {10^{ - 16}}$12.042.1667$ \times {10^{ - 29}}$16.01
      162.7591$ \times {10^{ - 13}}$6.011.4223$ \times {10^{ - 19}}$10.014.9470$ \times {10^{ - 19}}$12.013.2973$ \times {10^{ - 34}}$16.00
      324.3152$ \times {10^{ - 15}}$6.001.3874$ \times {10^{ - 22}}$10.001.2058$ \times {10^{ - 22}}$12.005.0279$ \times {10^{ - 39}}$16.00
      $\bar \rho \approx $6101216

      例2. 反复修正结点

      微分方程同例1。仅考虑单一线性元($m = 1$, ${N_e} = 1$),其两个结点记作0 (${x_0} = 0$, $u_{\rm{0}}^h = 0$)和1 (${x_{\rm{1}}} = {\rm{1}}$, $u_1^h$待求)。本例采用EEP简约格式,用第2.3节算法反复修正右端点位移$u_1^h$,逐步提升其精度;下面给出具体计算过程。

      对于初始有限元解,容易得到刚度系数${k_{{\rm{11}}}} = {{11} / 6}$和结点荷载${P_{\rm{1}}} = {1 / 2}$,则:

      $$ {u^h} = {{3x} / {11}},\;u_{\rm{1}}^h = {3 / {11}} $$ (21)

      由式(10a),得EEP解为:

      $$ {u^*} = {{x(13 - 8x + {x^2})} / {22}},\;u_{\rm{1}}^* = u_{\rm{1}}^h = {3 / {11}} $$ (22)

      可见,从${u^h}$${u^*}$,解函数从线性提升到三次函数。用残余荷载$1 - L\; {u^*}$代替原始荷载1,重新计算结点荷载得到$P_{\rm{1}}^* = {1 / {40}}$,可得:

      $$ \varepsilon _{\rm{1}}^h = {3 / {220}},\;\tilde u_{\rm{1}}^h = u_{\rm{1}}^h + \varepsilon _{\rm{1}}^h{{ = 63} / {220}} $$ (23)

      再用式(10a)计算EEP解得:

      $$ \begin{split} & {\varepsilon ^*} = - \frac{{23}}{{330}}x + \frac{{73}}{{440}}{x^2} - \frac{{29}}{{440}}{x^3} - \frac{5}{{264}}{x^4} + \frac{1}{{440}}{x^5},\\& \varepsilon _{\rm{1}}^* = \varepsilon _{\rm{1}}^h = {3 / {220}}\\[-16pt] \end{split} $$ (24)

      后续重复计算的过程及其结果见表3。可以看到,本法所估计的结点位移误差精度相当高(如估计值为$\varepsilon _1^h = 0.0136363636$,而真实值为$e_1^h = 0.0116050160$),且随着迭代次数增加,结点位移趋近于精确解。

      表 3  反复修正结点位移

      Table 3.  Repeated nodal value updating

      结点1有限元解修正的结点值结点误差
      $u_{\rm{1}}^h$ 0.2727272727 0.2727272727 0.0116050160
      $\varepsilon _{\rm{1}}^h$ 0.0136363636 0.2863636363 −0.0020313477
      $\delta \varepsilon _{\rm{1}}^h$ −0.0021467926 0.2842168437 0.0001154449
      ${\delta ^2}\varepsilon _{\rm{1}}^h$ 0.0000908271 0.2843076709 0.0000246178
      ${\delta ^{\rm{3}}}\varepsilon _{\rm{1}}^h$ 0.0000300643 0.2843377351 −0.0000054464
      ${\delta ^{\rm{4}}}\varepsilon _{\rm{1}}^h$ −0.0000055624 0.2843321727 0.0000001159
      ··· ··· ··· ···
      精确解 0.2843322887

      表4是4个线性元的结点误差,可以看到网格加密后精度随之提高。而且,尽管没有展示,单元次数更高也会得到更高的精度。

      例3. 初值问题的附加精度

      本例中,考虑简谐荷载作用下的有阻尼的单自由度体系的运动方程:

      $$ \left\{ \begin{aligned} & m\,u'' + c\; u' + \; \; k\,u = \sin (\omega \,t),\quad 0 < t < T \\& u(0) = 0,\quad \quad u'(0) = 1 \end{aligned} \right. $$ (25)

      式中:无量纲的系数为$m = k = 1$$c\; = 0.04$$\omega = 0.2$$T = 125$。此问题的精确解(如图1所示)为:

      $$ \begin{split} & u(t) = {e^{ - \zeta t}}\left( {A\cos ({\omega _D}t) + B\sin ({\omega _D}t)} \right) + \\ &\qquad\quad C\sin (\omega \,t) + D\cos (\omega \,t) ,\\ & C = \frac{{1 - {\omega ^2}}}{{{{(1 - {\omega ^2})}^2} + {{(2\zeta \omega )}^2}}},\\& D = \frac{{ - 2\zeta \omega }}{{{{(1 - {\omega ^2})}^2} + {{(2\zeta \omega )}^2}}},\\& {\omega _D} = \sqrt {1 - {\zeta ^2}} ,\;A = - D,\\& B = (1 - \omega \,C - \zeta D)/{\omega _D} \end{split} $$ (26)

      对于初值问题,没必要形成整体刚度矩阵,有限元解可以从左到右逐个单元积分求解,因此结点位移修正算法可以逐单元进行,而不用像边值问题待整体求解之后再作修正。如此,对于初值问题,便出现一个增强的修正方法,其特点是在计算单元右端点有限元位移解$u_2^h$时,其左端点的初始位移,可以立即采用对有限元解$u_1^h$修正后的结点值$\tilde u_1^h$。这个逐单元修正的方法,不仅更高效,而且结点精度得到进一步的显著提高。

      表 4  4个线性元的结点误差

      Table 4.  Errors of nodal values of four linear elements

      网格${N_e} = 4$,$m = 1$
      结点 $i = 1$ $i = 2$ $i = 3$ $i = 4$
      $\varepsilon _i^h$ −0.1802$ \times {10^{ - 3}}$ −0.2242$ \times {10^{ - 3}}$ 0.2363$ \times {10^{ - 4}}$ 0.6194$ \times {10^{ - 3}}$
      $\varepsilon _i^h + \delta \varepsilon _i^h$ −0.1776$ \times {10^{ - 3}}$ −0.2198$ \times {10^{ - 3}}$ 0.2017$ \times {10^{ - 4}}$ 0.6153$ \times {10^{ - 3}}$
      $e_i^h$ −0.1775$ \times {10^{ - 3}}$ −0.2197$ \times {10^{ - 3}}$ 0.2022$ \times {10^{ - 4}}$ 0.6146$ \times {10^{ - 3}}$

      图  1  例3的精确解

      Figure 1.  Exact solution of Example 3

      分别使用单元长度(时间步长)为$h = 0.4,\;0.2,\;0.1$的线性元计算。表5给出三个求解方案的最大结点位移误差,即:1)有限元解$u_i^h$;2)整体修正解$u_i^h + \varepsilon _i^h$(同边值问题);3)逐单元修正解${(u_i^h + \varepsilon _i^h)^{*}}$(仅限于初值问题)。可以看到,尽管方案2)和方案3)的收敛阶相同$O({h^4})$,但是方案3)的逐单元计算使得精度有翻倍提升:大约是方案2)的16倍,是方案1)的256倍。图2~图4分别给出三个方案的结点位移误差(固定单元步长$h = 0.2$),易见逐单元法具有显著的精度优势。

      图  2  有限元解误差

      Figure 2.  Errors of FE solution

      图  3  整体修正误差

      Figure 3.  Errors of globally updated solution

      图  4  逐单元修正误差

      Figure 4.  Errors of element-wise updated solution

      表 5  结点误差和收敛阶($m = 1$,简约格式)

      Table 5.  Nodal errors and convergence orders ($m = 1$, simplified form)

      单元步长$h$有限元解$u_i^h$整体修正解$u_i^h + \varepsilon _i^h$逐单元修正解${(u_i^h + \varepsilon _i^h)^*}$
      $\mathop {\max }\limits_i \left| {u_i^{} - u_i^h} \right|$$\rho $$\mathop {\max }\limits_i \left| {u_i^{} - u_i^h - \varepsilon _i^h} \right|$$\rho $$\mathop {\max }\limits_i \left| {u_i^{} - {{(u_i^h + \varepsilon _i^h)}^*}} \right|$$\rho $
      0.496.878$ \times {10^{ - 3}}$23.437$ \times {10^{ - 3}}$1.449$ \times {10^{ - 3}}$
      0.224.400$ \times {10^{ - 3}}$1.991.484$ \times {10^{ - 3}}$3.980.095$ \times {10^{ - 3}}$3.92
      0.16.102$ \times {10^{ - 3}}$2.000.093$ \times {10^{ - 3}}$4.000.006$ \times {10^{ - 3}}$3.98
      $\bar \rho \approx $244
    • 虽然本文的方法和算例主要针对二阶常微分方程,但是所提出的方法也适用于其它一维问题,如其它阶常微分方程问题以及常微分方程组问题。此外,所提出的方法有潜力通过所谓的逐维修正方法扩展到二维乃至三维有限元分析[10-11]

参考文献 (15)

目录

    /

    返回文章
    返回