留言板

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

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

多边形比例边界有限元非线性高效分析方法

李佳龙 李钢 余丁浩

李佳龙, 李钢, 余丁浩. 多边形比例边界有限元非线性高效分析方法[J]. 工程力学, 2020, 37(9): 8-17. doi: 10.6052/j.issn.1000-4750.2019.10.0634
引用本文: 李佳龙, 李钢, 余丁浩. 多边形比例边界有限元非线性高效分析方法[J]. 工程力学, 2020, 37(9): 8-17. doi: 10.6052/j.issn.1000-4750.2019.10.0634
Jia-long LI, Gang LI, Ding-hao YU. POLYGON SCALED BOUNDARY FINITE ELEMENT NONLINEAR EFFICIENT ANALYSIS METHOD[J]. Engineering Mechanics, 2020, 37(9): 8-17. doi: 10.6052/j.issn.1000-4750.2019.10.0634
Citation: Jia-long LI, Gang LI, Ding-hao YU. POLYGON SCALED BOUNDARY FINITE ELEMENT NONLINEAR EFFICIENT ANALYSIS METHOD[J]. Engineering Mechanics, 2020, 37(9): 8-17. doi: 10.6052/j.issn.1000-4750.2019.10.0634

多边形比例边界有限元非线性高效分析方法

doi: 10.6052/j.issn.1000-4750.2019.10.0634
基金项目: 国家重点研发计划项目(2018YFC1504303);国家自然科学基金项目(51878112);辽宁省“兴辽英才计划”项目(XLYC1902043);大连市高层次人才创新支持计划项目(2017RD04)
详细信息
    作者简介:

    李佳龙(1991−),男,四川遂宁人,博士生,主要从事结构非线性分析等研究(E-mail: wekalee@163.com)

    余丁浩(1989−),男,河北邯郸人,博士后,主要从事结构非线性分析等研究(E-mail: ydh@dlut.edu.cn)

    通讯作者: 李 钢(1979−),男,辽宁葫芦岛人,教授,博士,博导,主要从事结构工程抗震等研究(E-mail: gli@dlut.edu.cn)
  • 中图分类号: O346.1;O302

POLYGON SCALED BOUNDARY FINITE ELEMENT NONLINEAR EFFICIENT ANALYSIS METHOD

  • 摘要: 比例边界有限元法作为一种高精度的半解析数值求解方法,特别适合于求解无限域与应力奇异性等问题,多边形比例边界单元在模拟裂纹扩展过程、处理局部网格重剖分等方面相较于有限单元法具有明显优势。目前,比例边界有限元法更多关注的是线弹性问题的求解,而非线性比例边界单元的研究则处于起步阶段。该文将高效的隔离非线性有限元法用于比例边界单元的非线性分析,提出了一种高效的隔离非线性比例边界有限元法。该方法认为每个边界线单元覆盖的区域为相互独立的扇形子单元,其形函数以及应变-位移矩阵可通过半解析的弹性解获得;每个扇形区的非线性应变场通过设置非线性应变插值点来表达,引入非线性本构关系即可实现多边形比例边界单元高效非线性分析。多边形比例边界单元的刚度通过集成每个扇形子单元的刚度获取,扇形子单元的刚度可采用高斯积分方案进行求解,其精度保持不变。由于引入了较多的非线性应变插值点,舒尔补矩阵维数较大,该文采用Woodbury近似法对隔离非线性比例边界单元的控制方程进行求解。该方法对大规模非线性问题的计算具有较高的计算效率,数值算例验证了算法的正确性以及高效性,将该方法进行推广,对实际工程分析具有重要意义。
  • 图  1  多边形比例边界单元

    Figure  1.  Polygon scaled boundary element

    图  2  多边形单元高斯点布置方案

    Figure  2.  Gaussian points layout scheme of polygon element

    图  3  隔离非线性比例边界有限元法分析流程

    Figure  3.  Analysis flow of IS-SBFEM

    图  4  悬臂梁模型

    Figure  4.  Cantilever beam model

    图  5  多边形网格

    Figure  5.  Polygon mesh

    图  6  有限元网格

    Figure  6.  Finite element mesh

    图  7  竖向荷载-位移曲线

    Figure  7.  Curves of the vertical load vs. displacement

    图  8  计算自由度曲线

    Figure  8.  Calculate DOFs curves

    图  9  几何模型与多边形网格

    Figure  9.  Geometry model and polygon mesh

    图  10  Koyna地震动记录

    Figure  10.  Seismic acceleration records of Koyna

    图  11  坝顶横向位移反应

    Figure  11.  Horizontal displacement responses of dam crest

    图  12  坝顶竖向位移反应

    Figure  12.  Vertical displacement responses of dam crest

    图  13  坝顶横向速度反应

    Figure  13.  Horizontal velocity responses of dam crest

    图  14  坝顶竖向速度反应

    Figure  14.  Vertical velocity responses of dam crest

    图  15  坝顶横向加速度反应

    Figure  15.  Horizontal acceleration responses of dam crest

    图  16  坝顶横向加速度反应

    Figure  16.  Vertical acceleration responses of dam crest

    图  17  计算自由度曲线

    Figure  17.  Calculate DOFs curves

    表  1  扇形单元积分点坐标

    Table  1.   Integral point coordinates of sector element

    iξηω
    1${{{\rm{(}}{{1 - 1} / {\sqrt 3 }}{\rm{)}}} / 2}$${{ - 1} / {\sqrt 3 }}$1/2
    2${{{\rm{(}}{{1 + 1} / {\sqrt 3 }}{\rm{)}}} / 2}$${{ - 1} / {\sqrt 3 }}$1/2
    3${{{\rm{(}}{{1 + 1} / {\sqrt 3 }}{\rm{)}}} / 2}$${1 / {\sqrt 3 }}$1/2
    4${{{\rm{(}}{{1 - 1} / {\sqrt 3 }}{\rm{)}}} / 2}$${1 / {\sqrt 3 }}$1/2
    下载: 导出CSV

    表  2  A点处的位移

    Table  2.   Displacement of node A

    模型横向位移/mm相对误差/(%)竖向位移/mm相对误差/(%)
    FEM-a0.340312.7202.31762.255
    FEM-b0.38321.7202.56930.864
    FEM-c0.38992.5917
    半解析SBFEM0.38311.7202.56810.911
    数值SBFEM0.38321.7202.56900.876
    下载: 导出CSV

    表  3  控制方程计算时间

    Table  3.   Calculate time of the governing equation

    计算方法方程求解时间/s
    SBFEM 3.19
    近似IS-SBFEM 8.56
    精确IS-SBFEM202.96
    下载: 导出CSV

    表  4  控制方程计算时间

    Table  4.   Calculate time of governing equation

    计算方法方程求解时间/s
    SBFEM410.60
    近似IS-SBFEM 58.57
    精确IS-SBFEM
    下载: 导出CSV
  • [1] Wolf J P. The scaled boundary finite element method [M]. West Sussex: Wily, 2003.
    [2] Song C M. The scaled boundary finite element method: Introduction to theory and implementation [M]. New Jersey: John Wiley & Sons Ltd, 2018.
    [3] Deeks A J, Wolf J P. A virtual work derivation of the scaled boundary finite-element method for elastostatics [J]. Computational Mechanics, 2002, 28(6): 489 − 504.
    [4] Natarajan S, Ooi E T, Chiong I, et al. Convergence and accuracy of displacement based finite element formulations over arbitrary polygons: Laplace interpolants, strain smoothing and scaled boundary polygon formulation [J]. Finite Element in Analysis and Design, 2014, 85: 101 − 122. doi:  10.1016/j.finel.2014.03.006
    [5] Jia Y M, Li C J, Zhang Y, et al. The high-order completeness analysis of the scaled boundary finite element method [J]. Computer Methods in Applied Mechanics and Engineering, 2020, 362(4): 112867.
    [6] Bazyar M H, Song C M. A continued-fraction-based high-order transmitting boundary for wave propagation in unbounded domains of arbitrary geometry [J]. International Journal for Numerical Methods in Engineering, 2008, 74: 209 − 237. doi:  10.1002/nme.2147
    [7] 钟红, 暴艳利, 林皋. 基于多边形比例边界有限元的重力坝裂纹扩展过程模拟[J]. 水利学报, 2014, 45(增刊 1): 30 − 37.

    Zhong Hong, Bao Yanli, Lin Gao. Modeling of crack propagation of gravity dams based on scaled boundary polygons [J]. Journal of Hydraulic Engineering, 2014, 45(Suppl 1): 30 − 37. (in Chinese)
    [8] 李志远, 李建波, 林皋. 局部复杂场地Rayleigh波传播特性分析的SBFEM型研究[J]. 岩土力学, 2018, 39(11): 4242 − 4250.

    Li Zhiyuan, Li Jianbo, Lin Gao. Research on influence of partial terrain to scattering of Rayleigh wave based on SBFEM [J]. Rock and Soil Mechanics, 2018, 39(11): 4242 − 4250. (in Chinese)
    [9] Lin G, Li Z Y, Li J B. A substructure replacement technique for the numerical solution of wave scattering problem [J]. Soil Dynamics and Earthquake Engineering, 2018, 111: 87 − 97. doi:  10.1016/j.soildyn.2018.04.031
    [10] Ooi E T, Song C M, Tin-Loi F. A scaled boundary polygon formulation elasto-plastic analyses [J]. Computer Methods in Applied Mechanics and Engineering, 2014, 268: 905 − 937.
    [11] Chen K, Zou D G, Kong X J, et al. A novel nonlinear solution for the polygon scaled boundary finite element method and its application to geotechnical structures [J]. Computers and Geotechnics, 2017, 82: 201 − 210. doi:  10.1016/j.compgeo.2016.09.013
    [12] Li G, Yu D H. Efficient inelasticity-separated finite-element method for material nonlinearity analysis [J]. Journal of Engineering Mechanics, 2018, 144(4): 04018008. doi:  10.1061/(ASCE)EM.1943-7889.0001426
    [13] Akgün M A, Garcelon J H, Haftka R T. Fast exact linear and non-linear structural reanalysis and the Sherman–Morrison–Woodbury formulas [J]. International Journal for Numerical Methods in Engineering, 2001, 50(7): 1587 − 1606. doi:  10.1002/nme.87
    [14] Li G, Li J L, Yu L, et al. Improved Woodbury approximation approach for inelasticity-separated solid model analysis [J]. Soil Dynamics and Earthquake Engineering, 2020, 129: 105926. doi:  10.1016/j.soildyn.2019.105926
    [15] Li G, Jia S, Yu D H, et al. Woodbury approximation method for structural nonlinear analysis [J]. Journal of Engineering Mechanics, 2018, 144(7): 04018052. doi:  10.1061/(ASCE)EM.1943-7889.0001464
    [16] Li G, Yu D H, Li H N. Seismic response analysis of reinforced concrete frames using inelasticity-separated fiber beam-column model [J]. Earthquake Engineering & Structural Dynamics, 2018, 47(5): 1291 − 1308.
    [17] 李佳龙, 李钢, 于龙. 隔离非线性平面应变单元模型及其在Drucker-Prager模型中的应用[J]. 岩土力学, 2020, 41(5): 1492 − 1501.

    Li Jialong, Li Gang, Yu Long. The inelasticity-separated plane strain element model and its application in drucker-prager model [J]. Rock and Soil Mechanics, 2020, 41(5): 1492 − 1501. (in Chinese)
    [18] 李佳龙, 李钢, 李宏男. 基于隔离非线性的实体单元模型与计算效率分析[J]. 工程力学, 2019, 36(9): 40 − 49.

    Li Jialong, Li Gang, Li Hongnan. The inelastcity-separated solid element model and computational efficiency analysis [J]. Engineering Mechanics, 2019, 36(9): 40 − 49. (in Chinese)
    [19] 黄羽立, 陆新征, 叶列平, 等. 基于多点位移控制的推覆分析算法[J]. 工程力学, 2011, 28(2): 18 − 23.

    Huang Yuli, Lu Xinzheng, Ye Lieping, et al. A pushover analysis algorithm based on multiple point constraints [J]. Engineering Mechanics, 2011, 28(2): 18 − 23. (in Chinese)
    [20] Lee J, Fenves G L. A plastic-damage concrete model for earthquake analysis of dams [J]. Earthquake Engineering & Structural Dynamics, 1998, 27(9): 937 − 956.
    [21] 杨强, 杨晓君, 陈新. 基于D-P准则的理想弹塑性本构关系积分研究[J]. 工程力学, 2005, 22(4): 15 − 19. doi:  10.3969/j.issn.1000-4750.2005.04.004

    Yang Qiang, Yang Xiaojun, Chen Xin. On integration algorihms for perfect plasticity based on drucker-prager criterion [J]. Engineering Mechanics, 2005, 22(4): 15 − 19. (in Chinese) doi:  10.3969/j.issn.1000-4750.2005.04.004
  • [1] 苏璞, 李钢, 余丁浩.  基于子结构的Woodbury非线性分析方法 . 工程力学, 2020, 37(5): 26-35. doi: 10.6052/j.issn.1000-4750.2019.07.0419
    [2] 李钢, 吕志超, 余丁浩.  隔离非线性分层壳有限单元法 . 工程力学, 2020, 37(3): 18-27. doi: 10.6052/j.issn.1000-4750.2019.01.0189
    [3] 张微敬, 张晨骋.  钢筋套筒挤压连接的预制RC柱非线性有限元分析 . 工程力学, 2018, 35(S1): 67-72. doi: 10.6052/j.issn.1000-4750.2018.01.S006
    [4] 钟红, 林皋, 李红军.  坝基界面在非线性水压力驱动下的非线性断裂过程模拟 . 工程力学, 2017, 34(4): 42-48. doi: 10.6052/j.issn.1000-4750.2015.10.0817
    [5] 李上明, 吴连军.  基于连分式与有限元法的坝库耦合瞬态分析方法 . 工程力学, 2016, 33(4): 9-16. doi: 10.6052/j.issn.1000-4750.2014.09.0803
    [6] 陈政清.  梁杆结构几何非线性有限元的数值实现方法 . 工程力学, 2014, 31(6): 42-52. doi: 10.6052/j.issn.1000-4750.2013.05.ST08
    [7] 寇佳亮, 梁兴文, 钱 磊, 邓明科.  高强混凝土剪力墙非线性分析及抗震性能参数研究 . 工程力学, 2013, 30(2): 232-239. doi: 10.6052/j.issn.1000-4750.2011.06.0381
    [8] 李上明.  基于比例边界有限元法动态刚度矩阵的坝库耦合分析方法 . 工程力学, 2013, 30(2): 313-317. doi: 10.6052/j.issn.1000-4750.2011.08.0566
    [9] 李上明.  基于比例边界有限元法和有限元法的水下结构瞬态分析方法 . 工程力学, 2013, 30(11): 42-46. doi: 10.6052/j.issn.1000-4750.2012.07.0532
    [10] 朱明亮, 董石麟.  向量式有限元在索穹顶静力分析中的应用 . 工程力学, 2012, 29(8): 236-242. doi: 10.6052/j.issn.1000-4750.2010.12.0874
    [11] 秦 剑, 黄克服, 张清东.  几何非线性样条有限元法 . 工程力学, 2011, 28(增刊I): 1-004.
    [12] 于琦, 孟少平, 吴京, 郑开启.  预应力混凝土结构组合式非线性分析模型 . 工程力学, 2011, 28(11): 130-137.
    [13] 王德才, 叶献国, 连 星, 常 磊.  基于纤维模型叠合板式剪力墙非线性分析 . 工程力学, 2010, 27(增刊I): 164-167.
    [14] 桂国庆, 王玉娥.  铝合金单层球面网壳的非线性稳定分析 . 工程力学, 2006, 23(增刊Ⅱ): 0-035,.
    [15] 尚守平, 曾令宏, 彭晖.  钢丝网复合砂浆加固混凝土受弯构件非线性分析 . 工程力学, 2005, 22(3): 118-125.
    [16] 周奎, 方立新, 宋启根.  组合结构梁柱COMP单元非线性分析 . 工程力学, 2001, 18(4): 67-72.
    [17] 陈务军, 董石麟, 付功义, 周岱, 胡继军.  非线性分析广义增量法算法研究 . 工程力学, 2001, 18(3): 28-33,4.
    [18] 高向宇, 方鄂华.  组合墙结构非线性分析模型 . 工程力学, 1997, 14(2): 36-42.
    [19] 曲庆璋, 梁兴复.  弹性地基上自由矩形板的非线性动静态分析 . 工程力学, 1996, 13(3): 40-46.
    [20] 陈进, 陈祥福.  岩石工程渐进破坏分析的边界元法 . 工程力学, 1991, 8(4): 128-137.
  • 加载中
图(17) / 表 (4)
计量
  • 文章访问数:  68
  • HTML全文浏览量:  38
  • PDF下载量:  18
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-10-29
  • 修回日期:  2020-04-07
  • 网络出版日期:  2020-09-07
  • 刊出日期:  2020-09-25

多边形比例边界有限元非线性高效分析方法

doi: 10.6052/j.issn.1000-4750.2019.10.0634
    基金项目:  国家重点研发计划项目(2018YFC1504303);国家自然科学基金项目(51878112);辽宁省“兴辽英才计划”项目(XLYC1902043);大连市高层次人才创新支持计划项目(2017RD04)
    作者简介:

    李佳龙(1991−),男,四川遂宁人,博士生,主要从事结构非线性分析等研究(E-mail: wekalee@163.com)

    余丁浩(1989−),男,河北邯郸人,博士后,主要从事结构非线性分析等研究(E-mail: ydh@dlut.edu.cn)

    通讯作者: 李 钢(1979−),男,辽宁葫芦岛人,教授,博士,博导,主要从事结构工程抗震等研究(E-mail: gli@dlut.edu.cn)
  • 中图分类号: O346.1;O302

摘要: 比例边界有限元法作为一种高精度的半解析数值求解方法,特别适合于求解无限域与应力奇异性等问题,多边形比例边界单元在模拟裂纹扩展过程、处理局部网格重剖分等方面相较于有限单元法具有明显优势。目前,比例边界有限元法更多关注的是线弹性问题的求解,而非线性比例边界单元的研究则处于起步阶段。该文将高效的隔离非线性有限元法用于比例边界单元的非线性分析,提出了一种高效的隔离非线性比例边界有限元法。该方法认为每个边界线单元覆盖的区域为相互独立的扇形子单元,其形函数以及应变-位移矩阵可通过半解析的弹性解获得;每个扇形区的非线性应变场通过设置非线性应变插值点来表达,引入非线性本构关系即可实现多边形比例边界单元高效非线性分析。多边形比例边界单元的刚度通过集成每个扇形子单元的刚度获取,扇形子单元的刚度可采用高斯积分方案进行求解,其精度保持不变。由于引入了较多的非线性应变插值点,舒尔补矩阵维数较大,该文采用Woodbury近似法对隔离非线性比例边界单元的控制方程进行求解。该方法对大规模非线性问题的计算具有较高的计算效率,数值算例验证了算法的正确性以及高效性,将该方法进行推广,对实际工程分析具有重要意义。

English Abstract

李佳龙, 李钢, 余丁浩. 多边形比例边界有限元非线性高效分析方法[J]. 工程力学, 2020, 37(9): 8-17. doi: 10.6052/j.issn.1000-4750.2019.10.0634
引用本文: 李佳龙, 李钢, 余丁浩. 多边形比例边界有限元非线性高效分析方法[J]. 工程力学, 2020, 37(9): 8-17. doi: 10.6052/j.issn.1000-4750.2019.10.0634
Jia-long LI, Gang LI, Ding-hao YU. POLYGON SCALED BOUNDARY FINITE ELEMENT NONLINEAR EFFICIENT ANALYSIS METHOD[J]. Engineering Mechanics, 2020, 37(9): 8-17. doi: 10.6052/j.issn.1000-4750.2019.10.0634
Citation: Jia-long LI, Gang LI, Ding-hao YU. POLYGON SCALED BOUNDARY FINITE ELEMENT NONLINEAR EFFICIENT ANALYSIS METHOD[J]. Engineering Mechanics, 2020, 37(9): 8-17. doi: 10.6052/j.issn.1000-4750.2019.10.0634
  • 比例边界有限元法(Scaled boundary finite element method,简称SBFEM)是由Wolf和Song[1-2]提出的一种半解析的数值求解方法,最早用于解决土-结构相互作用问题中的无限域动力刚度的计算[3]。多边形比例边界有限元法[4]结合了有限单元法与边界元法的优势,具有较高的计算精度以及收敛性[5],目前已运用到众多工程问题的数值分析中[6-9]。多边形比例边界单元由于其网格的灵活性,在模拟裂纹扩展过程、处理局部网格重剖分等方面相较于有限元单元具有更明显的优势。目前,比例边界有限元法更多关注的是线弹性问题的求解,对于非线性问题的求解则研究较少。通过对多边形比例边界单元的弹性特性进行分析,可获得与有限单元法等价的形函数以及应变-位移矩阵,且对于线弹性问题与非线性问题均适用。基于此,Ooi等[10]开创了多边形比例边界单元的非线性分析,并将其运用于弹塑性带裂纹问题的求解,其中单元的弹塑性本构矩阵采用最小二乘法进行数值拟合,但该方法计算过程复杂且要求多边形单元尺寸足够小才能够精确表达塑性区的变化。Chen等[11]提出在每个线单元覆盖的扇形区域引入3个数值积分点来进行多边形比例边界单元的非线性分析,具有较高的计算精度。无论是Ooi等[10]亦或是Chen[11]等所提出的多边形比例边界有限元非线性分析过程,两者的基本思路均是:先计算多边形比例边界单元的弹塑性刚度矩阵,然后再集成整个求解区域的总体刚度矩阵。然而,同传统有限单元法一样,均是对总体刚度矩阵进行分解求解从而获得单元的位移响应,但大规模刚度矩阵的分解依旧是限制非线性问题高效求解最主要的因素,因而有必要提出一种针对多边形比例边界单元的高效求解方法。

    隔离非线性有限元法(Inelasticity-separated finite element method,简称IS-FEM)是由Li和Yu[12]基于有限单元法的基本框架提出的一种结构非线性分析新方法,采用Woodbury公式[13]直接求解控制方程,实现了局部非线性问题的高效求解。该方法将单元应变分解为线弹性应变与非线性应变两部分,非线性应变通过在单元中设置非线性应变插值点以插值的形式来表示,从而将结构的非线性部分隔离开来。对于非局部非线性问题,采用Woodbury近似法[14-15]对控制方程进行求解,将非线性问题的迭代求解过程转换为多次的常刚度矩阵回代以及稀疏矩阵与向量的乘积,对大规模非线性问题的求解具有明显的效率优势。目前已成功运用于钢筋混凝土框架结构[16]、二维[17]以及三维[18]等问题的非线性分析,显示出其较高的计算效率。

    本文将隔离非线性有限元法的思想用于比例边界有限元非线性分析,提出一种高效的隔离非线性比例边界有限元法(Inelasticity-separated scaled boundary finite element method,简称IS-SBFEM)。主要思想是:将多边形比例边界单元所有扇形区域独立考虑,在每个扇形单元内设置若干个非线性应变插值点,从而建立相互独立的非线性应变场。每个扇形区域的径向自然坐标以及环向自然坐标分别为ξ=[0,1]、η=[−1,1],因此本文针对多边形比例边界单元的扇形区域提出一种新的数值积分方案,即:采用4点高斯积分方案对线性的扇形子单元进行数值积分。对于更高阶的比例边界单元,选取更多的数值积分点进行数值积分即可获得足够的精度。该数值求解方案对于线弹性问题以及非线性问题均适用,且精度保持不变。由于隔离非线性比例边界单元相较于隔离非线性有限元单元引入了更多的非线性应变插值点,即便是局部非线性问题其舒尔补矩阵维数也往往较大,因此建议直接采用Woodbury近似法对控制方程进行求解,但该方法对较小规模的非线性问题计算效率偏低。鉴于此,对于较小规模的非线性问题,建议采用多边形比例边界有限元法直接求解;对于大规模的非线性问题,采用隔离非线性比例边界有限元法求解具有更高的计算效率。基于MATLAB平台分别编制了比例边界有限元法以及隔离非线性比例边界有限元法两套非线性分析程序,数值算例验证了本文所提出的隔离非线性比例边界有限元法的正确性与高效性。

    • 在求解任意的平面问题时,可采用满足比例中心要求的多边形比例边界单元进行离散。图1为一典型的多边形比例边界单元,连接比例中心O与外边界线单元端部结点可将其划分为多个扇形子单元。边界线单元可采用一维线单元进行离散,其自然坐标为−1≤η≤1,因而边界线单元上的任意一点坐标(xb(η), yb(η))可表示为:

      $$\left\{ {\begin{aligned} & {x_b}(\eta ) = {{N}}(\eta ){{{x}}_b} \\& {y_b}(\eta ) = {{N}}(\eta ){{{y}}_b} \end{aligned}} \right.$$ (1)

      式中:xb=[x1, x2, ···, xm]Tyb=[y1, y2, ···, ym]T为该边界线单元上m个结点的坐标;N(η)为一维Gauss-Lobatto-Lagrange形函数。

      $${{N}}(\eta ) = [ {\begin{array}{*{20}{c}} {{N_1}(\eta )}&{{N_2}(\eta )}& \cdots &{{N_m}(\eta )} \end{array}} ]$$ (2)

      图  1  多边形比例边界单元

      Figure 1.  Polygon scaled boundary element

      在比例中心沿边界的径向方向引入自然坐标ξ,其坐标范围为0≤ξ≤1,其中:ξ=0为比例中心位置,ξ=1为边界位置。因而该求解域内任意点坐标(x(ξ,η), y(ξ,η))可表示为:

      $$\left\{ {\begin{aligned} & x(\xi ,\eta ) = \xi \cdot {x_b}(\eta ) = \xi \cdot {{N}}(\eta ){{{x}}_b} \\& y(\xi ,\eta ) = \xi \cdot {y_b}(\eta ) = \xi \cdot {{N}}(\eta ){{{y}}_b} \end{aligned}} \right.$$ (3)

      比例边界有限元法假定单元位移在环向方向为数值解,而在径向方向则为一解析解。将其在环向与径向两个方向上解耦,可表示为:

      $$u(\xi ,\eta ) = {{{N}}_u}(\eta )u(\xi )$$ (4)

      式中:Nu(η)为边界线单元的插值形函数;u(ξ)为径向方向的位移函数,u(ξ)的求解可参考比例边界有限元法相关文献[1-210-11],此处不再赘述。通过对多边形比例边界单元的弹性特性进行分析,可获得每个扇形区与有限元单元等价的单元形函数Φ(ξ,η)以及应变-位移矩阵B(ξ,η)。两者仅与单元本身形状相关而与材料的性质无关,对线弹性问题以及非线性问题均适用,分别为:

      $${{\varPhi }}(\xi ,\eta ) = {{{N}}_u}(\eta ){{{\psi }}_u}{\xi ^{ - {{{S}}_n}}}{{\psi }}_u^{ - 1}\qquad\qquad\qquad\qquad\;$$ (5)
      $${{B}}(\xi ,\eta ) = \left[ {{{{B}}_1}(\eta ){{{\psi }}_u}[ - {{{S}}_n}] + {{{B}}_2}(\eta ){{{\psi }}_u}} \right]{\xi ^{ - {{{S}}_n} - I}}{{\psi }}_u^{ - 1}$$ (6)

      式中:ψu为位移模态对应的转换矩阵;Sn是由单元的Hamilton矩阵的负特征值组成的对角矩阵;B1(η)与B2(η)为单元由笛卡尔坐标到局部坐标的转换矩阵。此时,单元的应变场ε(ξ,η)为:

      $${{\varepsilon}} (\xi ,\eta ) = {{B}}(\xi ,\eta ){{{u}}_b}$$ (7)

      式中,ub为多边形比例边界单元的结点位移。

    • 比例边界有限元法作为一种基于位移的有限元求解方法,采用最小位能原理建立其非线性求解的控制方程。对于满足平衡条件的变形体结构有:

      $$ \begin{split} \int_\varOmega {\delta {{{\varepsilon }}^{\rm{T}}}\Delta {{\sigma }}(\xi ,\eta ){\rm{d}}\varOmega = }& \int_\varGamma {\delta {{u}}_b^{\rm{T}}{{{f}}_t}{\rm{d}}\varGamma } + \int_\varOmega {\delta {{u}}_b^{\rm{T}}{{{f}}_b}{\rm{d}}\varOmega } - \\& \int_\varOmega {\delta {{{\varepsilon }}^{\rm{T}}}{{\sigma }}(\xi ,\eta ){\rm{d}}\varOmega } \\[-15pt] \end{split}$$ (8)

      式中:ftfb分别为单元边界载荷与体积载荷;σ(ξ,η)为单元上一个荷载步收敛后的应力;δε为虚位移场δub对应的虚应变场,且有δε=B(ξ,η)δub。单元的应力增量Δσ(ξ,η)可通过单元的弹塑性矩阵Dep以及结点位移增量Δub表示为:

      $$\Delta {{\sigma}} (\xi ,\eta ) = {{{D}}_{\rm{{ep}}}}{{B}}(\xi ,\eta )\Delta {{{u}}_b}$$ (9)

      将式(5)、式(9)代入式(8)所示的虚功方程,即可得到多边形比例边界单元非线性求解的控制方程,即:

      $${{{k}}_{\rm{ep}}}\Delta {{{u}}_b} = \Delta {{f}}$$ (10)

      式中,kep、Δf分别为多边形比例边界单元的弹塑性刚度矩阵以及外荷载增量。

      $${{{k}}_{\rm{ep}}} = \int_\varOmega {{{B}}{{(\xi ,\eta )}^{\rm{T}}}{{{D}}_{\rm{ep}}}{{B}}(\xi ,\eta )} {\rm{d}}\varOmega \qquad\qquad $$ (11)
      $$ \begin{split} \Delta {{f}} = &\int_\varGamma {{{\varPhi }}{{(\xi ,\eta )}^{\rm{T}}}{{{f}}_t}{\rm{d}}\varGamma } + \int_\varOmega {{{\varPhi }}{{(\xi ,\eta )}^{\rm{T}}}{{{f}}_b}{\rm{d}}\varOmega } - \\& \int_\varOmega {{{{B}}}(\xi ,\eta )^{\rm{T}}{{\sigma }}(\xi ,\eta ){\rm{d}}\varOmega } \end{split} $$ (12)

      多边形比例边界单元的刚度矩阵kep与外荷载Δf的具体计算步骤为:首先,通过式(5)、式(6)获得多边形单元每一个扇形区的形函数Φ(ξ,η)以及应变-位移矩阵B(ξ,η),然后,通过式(11)、式(12)计算该扇形区对单元刚度矩阵以及外荷载的贡献,最后,再集成该多边形单元的刚度矩阵以及外荷载向量[11]。通过对所有的多边形单元进行计算并按自由度组装获得计算域的总体刚度矩阵以及外荷载向量,即可得到整个计算域的控制方程:

      $$\left( {\sum\limits_{i = 1}^{{N_e}} {{{{k}}_{\rm{ep}}}} } \right)\Delta {{X}} = {{{K}}_{\rm{ep}}}\Delta {{X}} = \Delta {{F}}$$ (13)

      式中:Kep即为整个计算域的弹塑性刚度矩阵;ΔX与ΔF分别为总的结点位移增量与外荷载向量;Ne为整个计算域多边形单元个数。采用上述方法,可将半解析的比例边界有限元法转换成与有限单元法一样的纯数值求解方法,从而实现了多边形比例边界单元的非线性分析。

      对于式(13)中的刚度矩阵Kep同样具有稀疏性,可采用LDLT分解法直接进行求解。然而,对于大规模问题的求解其计算效率偏低,因而有必要提出一种针对多边形比例边界单元的高效求解方法。

    • 隔离非线性有限元法(IS-FEM)[12]作为一种高效的材料非线性分析方法,已在二维、三维等有限元问题中展现其较高的计算效率,在此将其推广至多边形比例边界单元的非线性分析,提出一种高效的隔离非线性比例边界有限元法(IS-SBFEM)。对于一个多边形比例边界单元,采用增量形式,每个扇形区域内任意点处的应变增量Δε(ξ,η)可通过该扇形区的应变矩阵B(ξ,η)以及单元边界结点位移增量Δub表示为:

      $$\Delta {{\varepsilon}} (\xi ,\eta ) = {{B}}(\xi ,\eta )\Delta {{{u}}_b}$$ (14)

      基于弹塑性力学的基本理论,单元内任意点处的应变增量Δε(ξ,η)均可分解为线弹性的应变增量Δε′(ξ,η)以及非线性应变增量Δε″(ξ,η),即:

      $$\Delta {{\varepsilon}} (\xi ,\eta ) = \Delta {{\varepsilon}} '(\xi ,\eta ) + \Delta {{\varepsilon}} ''(\xi ,\eta )$$ (15)

      应力增量Δσ(ξ,η)可通过弹性矩阵De与线弹性应变Δε′(ξ,η)表示为:

      $$\Delta {{\sigma}} (\xi ,\eta ) = {{{D}}_e}\Delta {{\varepsilon}} '(\xi ,\eta )$$ (16)

      将每个扇形区域独立考虑,通过在每个扇形区内设置l个非线性应变插值点,以插值形式来表示扇形区的非线性应变场。

      $$ \begin{split} \Delta {{\varepsilon}} ''(\xi ,\eta ) =& {{{C}}_s}\Delta {{{{\varepsilon}} }_s''} = \left[ {{{{C}}_1}\;{{{C}}_2}\; \cdots \;{{{C}}_i}\; \cdots \;{{{C}}_l}} \right] \times \\& {[ {\Delta {{{\varepsilon}}}_1''^{\rm{T}}\;\Delta {{{\varepsilon}}}_2''^{\rm{T}}\; \cdots \;\Delta {{{\varepsilon}}}_i''^{\rm{T}}\; \cdots \;\Delta {{{\varepsilon}}}_l''^{\rm{T}}} ]^{\rm{T}}} \end{split} $$ (17)

      式中:Cs为多边形比例边界单元第s个扇形子单元的非线性应变插值矩阵;$\Delta {{{\varepsilon}} ''_s}$是由插值点处的非线性应变$\Delta {{{\varepsilon}} ''_i}$组成的非线性应变向量。将式(14)、(15)以及式(17)代入式(16),此时单元的应力增量可通过结点位移以及插值点处的非线性应变表示为:

      $$\Delta {{\sigma}} (\xi ,\eta ) = {{{D}}_e}{{B}}(\xi ,\eta )\Delta {{{u}}_b} - {{{D}}_e}{{{C}}_s}\Delta {{{\varepsilon}} ''_s}$$ (18)

      在一个非线性荷载步中,应力增量可由线性化的胡克定律表示为Δσ(ξ,η)=DepΔε(ξ,η),因而应力增量又可表示为:

      $$\Delta {{\sigma}} (\xi ,\eta ){\rm{ = (}}{{D}}_{\rm{ep}}^{ - 1} - {{D}}_e^{ - 1}{{\rm{)}}^{ - 1}}{{{C}}_s}\Delta {{{\varepsilon}} ''_s}$$ (19)
    • 比例边界有限元法作为一种基于位移的有限元分析方法,采用虚功原理可建立隔离非线性比例边界有限元法的控制方程。

      $$\int_\varOmega {\delta {{(\Delta {{\varepsilon }}(\xi ,\eta ))}^{\rm{T}}}\Delta {{\sigma }}(\xi ,\eta )} {\rm{d}}\varOmega = \delta {(\Delta {{{u}}_b})^{\rm{T}}}\Delta {{f}}$$ (20)

      式中,δε(ξ,η))为虚结点位移δub)对应的虚应变,即:δε(ξ,η))=B(ξ,η)δub)。其中,虚应变δε(ξ,η))同样可以分解为虚线弹性应变δε′(ξ,η))以及虚非线性应变δε″(ξ,η))两部分,即δε(ξ,η)) = δε′(ξ,η)) + δε″(ξ,η)),且同样满足式(15)所示的非线性应变插值关系:δε″(ξ,η))=Csδ($\Delta {\varepsilon ''_s}$)。将式(18)、式(19)两式所示的应力增量关系代入式(20),即可得到隔离非线性比例边界单元的控制方程:

      $$\left[ {\begin{array}{*{20}{c}} {{{{k}}_e}}&{{{k'}}} \\ {{{{{k'}}}^{\rm{T}}}}&{{{k''}}} \end{array}} \right]\left\{ \begin{array}{l} \;\;\;\Delta {{{u}}_b} \\ - \Delta {{{{\varepsilon}} _p''}} \\ \end{array} \right\} = \left\{ \begin{array}{l} \Delta {{f}} \\ \;\;{\bf{0}} \\ \end{array} \right\}$$ (21)

      式中:

      $$ {{{k}}_e} = \int_\varOmega {{{B}}{{(\xi ,\eta )}^{\rm{T}}}{{{D}}_e}{{B}}(\xi ,\eta ){\rm{d}}\varOmega } \qquad $$ (22)
      $$ {{k'}} = \int_\varOmega {{{B}}{{(\xi ,\eta )}^{\rm{T}}}{{{D}}_e}{{C}}{\rm{d}}\varOmega } \qquad\qquad\;\; $$ (23)
      $$ {{k''}} = \int_\varOmega {{{{C}}^{\rm{T}}}{{{D}}_e}{{({{{D}}_e} - {{{D}}_{\rm{ep}}})}^{ - 1}}{{{D}}_e}{{C}}{\rm{d}}} \varOmega $$ (24)

      式中:ke为多边形比例边界单元的初始刚度矩阵;k′为非线性应变与结点力之间的系数矩阵;k″为包含材料非线性信息的分块对角矩阵,在计算过程中实时变化。$\Delta {{{\varepsilon}} ''_p}$为整个多边形单元的非线性应变向量;矩阵C为整个多边形比例边界单元的非线性应变插值矩阵。

      $${{C}}{\rm{ = }}\left[ {{{{C}}_1}\;{{{C}}_2}\; \cdots \;{{{C}}_s}\; \cdots \;{{{C}}_{{n_s}}}} \right]$$ (25)
    • 一个多边形单元可拆分为ns个扇形区域,因此,式(22)所示的多边形比例边界单元的刚度矩阵ke的计算过程为:先单独对每个扇形区域Ωs进行积分运算,然后再对所有的扇形区域进行集成。每个扇形子单元积分区域Ωs的自然坐标均为ξ=[0,1]与η=[−1,1],且有dΩ=|J(η)|dηdξ。因此,初始弹性刚度矩阵ke可表示为:

      $$ \begin{split} {{{k}}_e} &\!=\! \sum\limits_{s = 1}^{{n_s}} {{{{k}}_{e,s}}} \!=\! \sum\limits_{s = 1}^{{n_s}} {\left( {\int_{{\varOmega _s}} {{{B}}{{(\xi ,\eta )}^{\rm{T}}}{{{D}}_e}{{B}}(\xi ,\eta ){\rm{d}}\varOmega } } \right)} \! = \\& \sum\limits_{s = 1}^{{n_s}} {\left( \!{\int_0^1\! {\int_{ - 1}^1 \!{h{{B}}{{(\xi ,\eta )}^{\rm{T}}}{{{D}}_e}{{B}}(\xi ,\eta )\xi \left| {{{J}}(\eta )} \right|} } {\rm{d}}\eta {\rm{d}}\xi } \!\right)} \end{split}$$ (26)

      式中:ke,s为第s个扇形区域对整个多边形单元刚度的贡献;h为单元厚度;|J(η)|为边界线单元的雅克比矩阵。在扇形区域内引入g个高斯积分点,矩阵ke,s可采用高斯积分方案进行数值积分,有:

      $${{{k}}_{e,s}} = \sum\limits_{i = 1}^g {h{\omega _i}{{B}}{{({\xi _i},{\eta _i})}^{\rm{T}}}{{{D}}_e}{{B}}({\xi _i},{\eta _i}){\xi _i}\left| {{{J}}({\eta _i})} \right|} $$ (27)

      式中,ωi为积分权重。当多边形比例边界单元的阶数较低时,可采用表1所示4点高斯积分方案进行计算,图2所示为多边形比例边界单元的数值积分点布置图,其中线单元上的高斯积分点仅用来计算多边形单元的形函数Φ(ξ,η)以及应变-位移矩阵B(ξ,η)。对于高阶的多边形比例边界单元,可设置更多的高斯积分点。

      该积分方案对式(11)的计算同样适用,仅需将式(26)中的弹性矩阵De更换为弹塑性矩阵Dep,即可得到多边形比例边界单元的弹塑性刚度矩阵kep。多边形比例边界单元的弹塑性矩阵Dep的计算可直接调用有限元程序中的弹塑性矩阵计算模块,从而完成多边形比例边界单元的非线性分析。

      表 1  扇形单元积分点坐标

      Table 1.  Integral point coordinates of sector element

      iξηω
      1${{{\rm{(}}{{1 - 1} / {\sqrt 3 }}{\rm{)}}} / 2}$${{ - 1} / {\sqrt 3 }}$1/2
      2${{{\rm{(}}{{1 + 1} / {\sqrt 3 }}{\rm{)}}} / 2}$${{ - 1} / {\sqrt 3 }}$1/2
      3${{{\rm{(}}{{1 + 1} / {\sqrt 3 }}{\rm{)}}} / 2}$${1 / {\sqrt 3 }}$1/2
      4${{{\rm{(}}{{1 - 1} / {\sqrt 3 }}{\rm{)}}} / 2}$${1 / {\sqrt 3 }}$1/2

      图  2  多边形单元高斯点布置方案

      Figure 2.  Gaussian points layout scheme of polygon element

      式(23)与式(24)所示的单元矩阵k′以及k″的计算同式(22)一致,均是先通过高斯积分方案计算每个扇形区域的值,然后再进行集成。在对第s个扇形区进行数值积分时,由于假定每个扇形区域非线性应变场相互独立,则此时仅有该扇形区插值矩阵Cs存在,其余扇形区域的插值矩阵均为零矩阵。因而单元矩阵k′与k″是以扇形区为单元的分块矩阵以及分块对角矩阵。

      $${{k'}} = [ {\begin{array}{*{20}{c|c|c|c|c|c}} {{{{{k}}}_1'}} & {{{{{k}}}_2'}} & \cdots & {{{{{k}}}_s'}} & \cdots & {{{{{k}}}_{{n_s}}'}} \end{array}} ]$$ (28)
      $${{k''}} = {\rm{diag}}\left( {{{{{k}}}_1''},{{{{k}}}_2''}, \cdots {{{{k}}}_s''} \cdots {{{{k}}}_{ns}''}} \right)\qquad$$ (29)

      每个扇形区域对应的矩阵${{{k}}_s'}$以及${{{k}}_s''}$分别为:

      $$ \begin{split} {{{{k}}}_s'} = & \int_{{\varOmega _s}} {{{B}}{{(\xi ,\eta )}^{\rm{T}}}{{{D}}_e}{{{C}}_s}{\rm{d}}\varOmega } = \\& \int_0^1 {\int_{ - 1}^1 {h{{B}}{{(\xi ,\eta )}^{\rm{T}}}{{{D}}_e}{{{C}}_s}\xi \left| {{{J}}(\eta )} \right|} } {\rm{d}}\eta {\rm{d}}\xi \end{split}\qquad\quad $$ (30)
      $$ \begin{split} {{{{k}}}_s''} =& \int_{{\varOmega _s}} {{{{C}}_s^{\rm{T}}}{{{D}}_e}{{({{{D}}_e} - {{{D}}_{\rm {ep}}})}^{ - 1}}{{{D}}_e}{{{C}}_s}{\rm{d}}\varOmega } = \\& \int_0^1 {\int_{ - 1}^1 {h{{{C}}_s^{\rm{T}}}{{{D}}_e}{{({{{D}}_e} - {{{D}}_{\rm {ep}}})}^{ - 1}}{{{D}}_e}{{{C}}_s}\xi \left| {{{J}}(\eta )} \right|{\rm{d}}\eta {\rm{d}}\xi } } \end{split} $$ (31)

      对于较低阶的比例边界单元,每个扇形区可采用图2所示的4个高斯积分点作为单元的非线性应变插值点。由于插值矩阵Cs为Kronecker delta函数,因此式(30)、式(31)中的矩阵分别是以插值点为单位的分块矩阵以及分块对角矩阵。

      $${{{k}}_s'} = [ {\begin{array}{*{20}{c|c|c|c}} {{{{{k}}}_{s1}'}} & {{{{{k}}}_{s2}'}} & {{{{{k}}}_{s3}'}} & {{{{{k}}}_{s4}'}} \end{array}} ]$$ (32)
      $${{{k}}_s''} = {\rm{diag}}({{{k}}_{s1}''},{{{k}}_{s2}''},{{{k}}_{s3}''},{{{k}}_{s4}''})$$ (33)

      其中,子矩阵分别为:

      $$\left\{ \!{\begin{aligned} & {{{{k}}}_{si}'} \!=\! h{\omega _i}{{B}}{({\xi _i},{\eta _i})^{\rm{T}}}{{{D}}_e}{\xi _i}\left| {{{J}}({\eta _i})} \right|,\;\;\;\;\;\;\;i\! =\! 1,2,3,4\\& {{{{k}}}_{si}''} \!= \!h{\omega _i}{{B}}{({\xi _i},{\eta _i})^{\rm{T}}}{{{D}}_e}{({{{D}}_e} - {{D}}_{\rm{ep}}^i)^{ - 1}}{{{D}}_e}{\xi _i}\left| {{{J}}({\eta _i})} \right| \end{aligned}} \right.$$ (34)

      式中,${{D}}_{\rm{ep}}^i $即为当前扇形区域第i个插值点处的弹塑性矩阵。

      将式(21)中非线性应变$\Delta {{{\varepsilon}} ''_p}$进行静力凝聚,可得:

      $$( {{{{k}}_e} - {{k'}}{{{{k''}}}^{ - 1}}{{{{k'}}}^{\rm{T}}}} )\Delta {{{u}}_b} = \Delta {{f}}$$ (35)

      可以证明,式(35)与式(10)两者完全等价,即矩阵kep=kekk−1kT,可以认为矩阵−kk−1kT是单元采用弹性刚度矩阵ke的修正矩阵。对整个计算域进行组装,可得整体控制方程为:

      $$\left[ {\begin{array}{*{20}{c}} {{{{K}}_e}}&{{{K'}}} \\ {{{{{K'}}}^{\rm{T}}}}&{{{K''}}} \end{array}} \right]\left\{ \begin{array}{l} \;\;\;\; \Delta {{X}} \\ - \Delta {{{{E}}}_p''} \\ \end{array} \right\} = \left\{ \begin{array}{l} \Delta {{F}} \\ \;\;{\bf{0}} \\ \end{array} \right\}$$ (36)

      将全局非线性应变$\Delta {{{E}}''_p}$进行静力凝聚,有:

      $$( {{{{K}}_e} - {{K'}}{{{{K''}}}^{ - 1}}{{{{K'}}}^{\rm{T}}}} )\Delta {{X}} = \Delta {{F}}$$ (37)

      同样的,可以证明其与式(13)完全等价。

    • 非线性分析时,对每个扇形区域的非线性应变插值点进行单元状态判定,并更新矩阵${{{K}}_m'}$${{{K}}_m''}$即可得到该荷载步的控制方程:

      $$\left[ {\begin{array}{*{20}{c}} {{{{K}}_e}}&{{{{{K}}}_m'}} \\ {{{K}}_m'^{\rm{T}}}&{{{{{K}}}_m''}} \end{array}} \right]\left\{ \begin{array}{l} \;\;\;\;\Delta {{{X}}_m} \\ - \Delta {{{{E}}}_{p,m}''} \\ \end{array} \right\} = \left\{ \begin{array}{l} \Delta {{{F}}_m} \\ \;\;\;{\bf{0}} \\ \end{array} \right\}$$ (38)

      基于Sherman-Morrison-Woodbury[13]公式对式(38)进行求解,有:

      $$\Delta {{{X}}_m} = ( {{{K}}_e^{ - 1} + {{K}}_e^{ - 1}{{{{K}}}_m'}{{K}}_{\rm{S}}^{ - 1}{{K}}_m'^{\rm{T}}{{K}}_e^{ - 1}} )\Delta {{{F}}_m}$$ (39)

      式中,KS为舒尔补矩阵,为一满阵。

      $${{{K}}_{\rm{S}}} = {{{K}}_m''} - {{K}}_m'^{\rm{T}}{{{K''}}^{ - 1}}{{{K}}_m'}$$ (40)

      舒尔补矩阵的维数即为隔离非线性法的计算自由度数,且也是决定该方法计算效率最主要的因素。采用Newton-Raphson方法进行非线性迭代求解,具有快速收敛特性。隔离非线性比例边界有限元法在迭代求解时,整个计算域的内力ΔFint可直接通过当前迭代步的结点位移增量ΔXm以及非线性应变向量$\Delta {{{E}}_{p,m}''}$获得,无须再对每个单元进行计算并集成。

      $$\Delta {{{F}}_{\rm{int}}} = {{{K}}_e}\Delta {{{X}}_m} - {{{K}}_m'}\Delta {{{E}}_{p,m}''}$$ (41)

      对于一个ns边形的隔离非线性比例边界单元,其内部共有4ns个非线性应变插值点,因而其舒尔补矩阵KS规模较大,直接求解的计算量可能比传统的刚度矩阵直接分解求解的计算量还大。在此,可采用Woodbury近似法[14-15](Woodbury approximation approach)对式(39)进行近似求解,整个求解过程仅为多次的常刚度矩阵回代以及稀疏矩阵与向量的乘积。该方法对于大规模非线性问题的求解具有明显优势,且规模越大其计算效率相对于直接的刚度矩阵分解求解算法效率更高。文献[17]已从理论上说明了对于平面有限元问题,当结点自由度数目超过2000时即可体现出近似隔离非线性有限元法的高效性。

      基于本文所提出的隔离非线性比例边界有限元模型,编制了相应的计算程序,整体分析流程如图3所示。

      图  3  隔离非线性比例边界有限元法分析流程

      Figure 3.  Analysis flow of IS-SBFEM

    • 图4为一悬臂梁模型,其几何尺寸如图所示,厚度为0.1 m。材料参数为:弹性模量E=2.0×105 MPa,泊松比为ν=0.2。采用多边形网格以及线性四边形有限元单元进行网格剖分,其中多边形网格共计109个单元以及264个结点,如图5所示。有限元网格采用三种网格尺寸,分别为80个单元与105个结点、320个单元与369个结点以及1280个单元与1377个结点,分别如图6(a)图6(b)以及图6(c)所示。

      图  4  悬臂梁模型

      Figure 4.  Cantilever beam model

      图  5  多边形网格

      Figure 5.  Polygon mesh

      图  6  有限元网格

      Figure 6.  Finite element mesh

    • 在自由端结点A处施加竖向荷载F=100 kN,并假定网格最密的有限元模型c为参考解。表2给出了结点A在不同网格模型下横向和竖向位移,可以看出采用本文提出的数值积分方案所得的计算结果与半解析的SBFEM结果基本一致,且与有限元模型b结果也基本一致。

      表 2  A点处的位移

      Table 2.  Displacement of node A

      模型横向位移/mm相对误差/(%)竖向位移/mm相对误差/(%)
      FEM-a0.340312.7202.31762.255
      FEM-b0.38321.7202.56930.864
      FEM-c0.38992.5917
      半解析SBFEM0.38311.7202.56810.911
      数值SBFEM0.38321.7202.56900.876
    • 假定该悬臂梁为双线性随动强化弹塑性模型,其初始屈服应力σy=360 MPa,切线模量Eh=0.002E。自由端结点A处施加往复竖向位移u,每个荷载步附加位移约束增量为d0=0.2 mm,共计2400个荷载步[19]。采用SBFEM以及IS-SBFEM两种求解方案,其中隔离非线性方案采用直接分解舒尔补矩阵的精确求解以及近似求解两种方法,共计四种求解方案。

      通过弹性分析可知,有限元模型b的刚度与多边形模型基本一致,因而选择模型b作为对比模型来验证非线性分析的精度。图7为两种网格模型结点A处的竖向荷载F与结点位移u的变化曲线,可以看出四种求解方案结果基本吻合,多边形单元的三种求解方案之间相对误差基本可以忽略。多边形模型与有限元模型计算结果一致,表明本文针对多边形比例边界单元所提出的数值积分方案以及IS-SBFEM的正确性。图8给出了IS-SBFEM与SBFEM两种求解方案的计算自由度曲线,比例边界有限元模型的刚度矩阵维数为510,而隔离非线性比例边界有限元模型的计算自由度随着非线性程度的增加与减少呈现高低变化。

      图  7  竖向荷载-位移曲线

      Figure 7.  Curves of the vertical load vs. displacement

      图  8  计算自由度曲线

      Figure 8.  Calculate DOFs curves

      表3给出了多边形网格三种求解方案控制方程的具体求解时间,可以看出SBFEM求解效率最高,近似IS-SBFEM法次之,精确IS-SBFEM最低。这是因为当问题规模较小时,刚度矩阵的分解与回代耗时两者差别不大,且仅需对刚度矩阵进行一次分解与回代即可得到当前迭代步的计算结果。近似的IS-SBFEM由于要进行多次的常刚度矩阵回代以及稀疏矩阵与向量的乘积,因而其计算效率相对于直接的刚度矩阵分解法反而更低。采用精确的IS-SBFEM进行求解时,由于舒尔补矩阵为高阶满阵且远大于刚度矩阵维数,因而其求解效率最低。

      表 3  控制方程计算时间

      Table 3.  Calculate time of the governing equation

      计算方法方程求解时间/s
      SBFEM 3.19
      近似IS-SBFEM 8.56
      精确IS-SBFEM202.96
    • 印度Koyna重力坝作为少数几个在强震中破坏且较有完整记录的重力坝之一,长期以来均被科研工作者奉为经典算例[20]。该坝长850 m,高103 m,符合平面应变模型几何特征,其具体几何尺寸如图9(a)所示。根据实际损伤以及相关文献可知,坝踵以及坡折面处均发生破坏,因而对这两处的网格进行局部加密,多边形有限元模型如图9(b)所示。其中坝踵以及坡折面处的网格细化图分别如图9(c)以及图9(d)所示,整个模型共计36904个多边形单元以及18753个结点。

      图  9  几何模型与多边形网格

      Figure 9.  Geometry model and polygon mesh

      坝体作用于基岩上,因而底部采用刚性边界。坝体材料参数分别选取为:弹性模量E=30.0 GPa,泊松比ρ=2630 kg/m3,混凝土抗拉强度为2.9 MPa,抗压强度为24.1 MPa。为验证本文所提算法的正确性以及高效性,选取Drucker-prager模型[21]模拟坝体的非线性效应。地震动激励采用Koyna地震波在横向和竖向两个方向加载,将地震动强度参数调整为0.7 g,时间间隔0.01 s,共计10 s,分别如图10(a)以及图10(b)所示。仅考虑坝体在自重以及地震荷载作用下反应,不考虑水压力对坝体反应的影响,结构阻尼比取0.05。

      图  10  Koyna地震动记录

      Figure 10.  Seismic acceleration records of Koyna

      图11~图16分别为坝顶横向与竖向位移、速度以及加速度时程对比曲线。其中有限元模型采用线性四边形等参单元进行网格划分,整个计算区域共剖分103081个结点以及102400个单元。可以看出,FEM、SBFEM以及IS-SBFEM三种算法的计算结果完全吻合,再次验证了本文所提出的隔离非线性比例边界有限元法的正确性。

      图  11  坝顶横向位移反应

      Figure 11.  Horizontal displacement responses of dam crest

      图  12  坝顶竖向位移反应

      Figure 12.  Vertical displacement responses of dam crest

      图  13  坝顶横向速度反应

      Figure 13.  Horizontal velocity responses of dam crest

      图  14  坝顶竖向速度反应

      Figure 14.  Vertical velocity responses of dam crest

      图17给出了该重力坝模型计算自由度曲线,虽然部分荷载步的计算自由度大于刚度矩阵维数,但大多数荷载步的计算自由度却很小甚至模型仍处于弹性状态。

      图  15  坝顶横向加速度反应

      Figure 15.  Horizontal acceleration responses of dam crest

      图  16  坝顶横向加速度反应

      Figure 16.  Vertical acceleration responses of dam crest

      图  17  计算自由度曲线

      Figure 17.  Calculate DOFs curves

      表4给出了两种比例边界有限元非线性分析方法控制方程的求解时间,可以看出SBFEM的计算时间远大于近似IS-SBFEM的计算时间,约为其7倍。该多边形模型其刚度矩阵维数为75406,其一次分解与回代的计算时间分别约为1.8 s与0.03 s。由于近似IS-SBFEM的一次迭代过程仅为多次常刚度矩阵回代以及稀疏矩阵与向量的乘积,对于大规模非线性问题,采用IS-SBFEM计算优势明显。由于存在高阶的舒尔补矩阵,精确的IS-SBFEM不再适用。可见,对于大规模问题,IS-SBFEM的计算效率远远大于SBFEM,此时IS-SBFEM应作为非线性分析的首选,既可在保证非线性计算精度的同时也具有极高的计算效率。

      表 4  控制方程计算时间

      Table 4.  Calculate time of governing equation

      计算方法方程求解时间/s
      SBFEM410.60
      近似IS-SBFEM 58.57
      精确IS-SBFEM
    • 基于多边形比例边界有限元法基本理论,将隔离非线性思想用于多边形比例边界有限元非线性分析。在每个边界线单元覆盖的扇形区域内设置若干个非线性应变插值点,并以插值的形式建立其非线性应变场,从而构造了隔离非线性多边形比例边界单元。采用虚功原理建立了隔离非线性比例边界单元的控制方程,单元的弹塑性矩阵可直接调用有限元程序的本构模块获得,从而实现了多边形比例边界有限元高效非线性分析。得到如下结论:

      (1)每个扇形区域其径向自然坐标以及环向自然坐标分别为ξ=[0,1]以及η=[−1,1],可采用标准高斯积分方案进行数值积分,且对于高阶多边形比例边界单元也同样适用。

      (2)多边形单元中引入较多的非线性应变插值点使得舒尔补矩阵维数较大,采用Woodbury近似法联合算法对控制方程进行求解,可以保证隔离非线性比例边界有限元法的计算精度与高效性。

      (3)对于大规模问题,建议采用隔离非线性比例边界有限元法进行分析,以便获得更高的计算效率。将该方法进行推广,对工程结构的非线性分析具有重要意义。

参考文献 (21)

目录

    /

    返回文章
    返回