留言板

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

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

基于共旋法与稳定函数的几何非线性平面梁单元

邓继华 谭建平 谭平 田仲初

邓继华, 谭建平, 谭平, 田仲初. 基于共旋法与稳定函数的几何非线性平面梁单元[J]. 工程力学, 2020, 37(11): 28-35. doi: 10.6052/j.issn.1000-4750.2020.01.0012
引用本文: 邓继华, 谭建平, 谭平, 田仲初. 基于共旋法与稳定函数的几何非线性平面梁单元[J]. 工程力学, 2020, 37(11): 28-35. doi: 10.6052/j.issn.1000-4750.2020.01.0012
Ji-hua DENG, Jian-ping TAN, Ping TAN, Zhong-chu TIAN. A GEOMETRIC NONLINEAR PLANE BEAM ELEMENT BASED ON COROTATIONAL FORMULATION AND ON STABILITY FUNCTIONS[J]. Engineering Mechanics, 2020, 37(11): 28-35. doi: 10.6052/j.issn.1000-4750.2020.01.0012
Citation: Ji-hua DENG, Jian-ping TAN, Ping TAN, Zhong-chu TIAN. A GEOMETRIC NONLINEAR PLANE BEAM ELEMENT BASED ON COROTATIONAL FORMULATION AND ON STABILITY FUNCTIONS[J]. Engineering Mechanics, 2020, 37(11): 28-35. doi: 10.6052/j.issn.1000-4750.2020.01.0012

基于共旋法与稳定函数的几何非线性平面梁单元

doi: 10.6052/j.issn.1000-4750.2020.01.0012
基金项目: 广东省自然科学基金项目(2015A030310141);中国博士后基金项目(2014M562154);国家自然科学基金项目(51478049);湖南省科技重大专项(2015GK1001-1);长沙理工大学土木工程优势特色重点学科创新性项目(2015ZDXK03);长沙理工大学研究生“实践创新与创业能力提升计划”项目(SJCX202001);长沙理工大学青年教师成长计划项目(2019QJCZ059)
详细信息
    作者简介:

    谭建平(1996−),男,湖南人,硕士生,主要从事大跨径桥梁施工监控研究(E-mail: 834115536@qq.com)

    谭 平(1973−),男,湖南人,研究员,博士,博导,主要从事结构抗震和防灾减灾研究(E-mail: ptan@gzu.edu.cn)

    田仲初(1963−),男,湖南人,教授,博士,博导,主要从事大跨径桥梁施工监控研究(E-mail: 382525361@qq.com)

    通讯作者: 邓继华(1975−),男,湖南冷水江人,副教授,博士,硕导,主要从事桥梁与结构非线性研究(E-mail: jihuadeng@sina.com)

A GEOMETRIC NONLINEAR PLANE BEAM ELEMENT BASED ON COROTATIONAL FORMULATION AND ON STABILITY FUNCTIONS

  • 摘要: 建立一个准确、高效的几何非线性梁单元对于描述杆系结构的非线性行为至关重要。该文基于共旋坐标法和稳定函数提出了一种几何非线性平面梁单元。该单元在形成中把变形和刚体位移分开,局部坐标系内采用稳定函数以考虑单元P-δ效应的影响,从局部坐标系到结构坐标系的转换则采用共旋坐标法以及微分以考虑几何非线性,给出了几何非线性平面梁单元在结构坐标系下的全量平衡方程和切线刚度矩阵;在此基础上根据带铰梁端弯矩为零的受力特征,导出了能考虑梁端带铰的单元切线刚度矩阵表达式。通过多个典型算例验证了算法与程序的正确性、计算精度和效率。
  • 图  1  变形前后平面梁单元

    Figure  1.  Plane beam element before and after deformation

    图  2  梁单元及微段

    Figure  2.  beam element and micro-segment

    图  3  投影增量法示意

    Figure  3.  Projection increment method

    图  4  Lee’s框架及荷载 /cm

    Figure  4.  Lee’s frame geometry and loading

    图  5  Lee’s框架

    Figure  5.  Lee’s frame

    图  6  端部承受集中荷载的悬臂梁

    Figure  6.  Cantilever subjected to concentrated load at free end

    图  7  对角点受拉力作用时的铰接方棱形框架

    Figure  7.  Diamond-shaped frame under a pare of opposite concentrated tensions in opposite angles

    表  1  悬臂梁自由端的位移

    Table  1.   Displacement of cantilever at free end

    荷载步数方法1个单元2个单元
    $u$$w$$\theta $$u$$w$$\theta $
    3152.33587.9181.45053.89383.4981.435
    256.76184.7781.478
    4152.33587.9181.45053.89383.4981.435
    2136.0593.2821.803
    5152.33587.9181.45053.89383.4981.435
    256.76184.7781.478
    6152.33587.9181.45053.89383.4981.435
    263.94493.2821.803
    7152.33587.9181.45053.89383.4981.435
    263.94493.2821.803
    解析解:$u$=55.5 $w$=81.06 $\theta $=1.430
    注:短横线表示该单元划分及荷载步条件下非线性计算失败。
    下载: 导出CSV

    表  2  计算误差

    Table  2.   Calculation error /(%)

    方法1个单元2个单元
    $u$$w$$\theta $$u$$w$$\theta $
    15.708.461.382.903.010.35
    215.2115.0826.052.274.593.31
    下载: 导出CSV

    表  3  铰接方棱形框架位移

    Table  3.   Displacement of diamond-shaped frame

    单元数$u/L$$w/L$${\theta _0}$最大误差/(%)
    10.46180.26247.6
    20.46190.25071.39976.9
    30.46520.24761.46382.6
    40.46510.24591.48161.5
    解析解0.46600.24381.5035
    注:在每根杆件划分成1个单元时,${\theta _0}$=0显然失真。
    下载: 导出CSV
  • [1] Yu J, Luo L, Li Y. Numerical study of progressive collapse resistance of RC beam-slab substructures under perimeter column removal scenarios [J]. Engineering Structures, 2018, 159: 14 − 27. doi:  10.1016/j.engstruct.2017.12.038
    [2] Qian Kai, Li Bing, Ma Jiaxing. Load-carrying mechanism to resist progressive collapse of RC buildings [J]. Journal of Structural Engineering, 2014, 141(2): 4014107 − 4014101.
    [3] 李钢, 吕志超, 余丁浩. 隔离非线性分层壳有限单元法[J]. 工程力学, 2020, 37(3): 18 − 27.

    Li Gang, Lv Zhichao, Yu Dinghao. The finite element model for inelasticity-separated multi-layer shell [J]. Engineering Mechanics, 2020, 37(3): 18 − 27. (in Chinese)
    [4] 李钢, 靳永强, 董志骞, 李宏男. 基于混合近似法的纤维梁单元非线性求解方法[J]. 土木工程学报, 2019, 52(6): 81 − 91.

    Li Gang, Jin Yongqiang, Dong Zhiqian, Li Hongnan. Nonlinear solution method for fiber beam element based on the hybrid approximations method [J]. China Civil Engineering Journal, 2019, 52(6): 81 − 91. (in Chinese)
    [5] 李佳龙, 李钢, 李宏男. 基于隔离非线性的实体单元模型与计算效率分析[J]. 工程力学, 2019, 36(9): 40 − 49.

    Li Jialong, Li Gang, Li Hongnan. The inelasticity-separated solid element model and computational efficiency analysis [J]. Engineering Mechanics, 2019, 36(9): 40 − 49. (in Chinese)
    [6] 王震, 李国强. 一种新型高效的几何非线性梁-柱单元[J]. 建筑钢结构进展, 2018, 20(3): 26 − 32.

    Wang Zhen, Li Guoqiang. A new efficient geometric nonlinear beam-column element [J]. Progress in Steel Building Structures, 2018, 20(3): 26 − 32. (in Chinese)
    [7] Rankin C C, Brogan F A. An element independent co-rotational procedure for the treatment of large rotation [J]. ASME J Pressure Vessel Technol, 1986, 108: 165 − 174. doi:  10.1115/1.3264765
    [8] Rankin C C, Nour-Omid B. The use of projectors to improve finite element performance [J]. Computers & Structures, 1988, 30: 257 − 267.
    [9] 史加贝, 刘铸永, 洪嘉振. 柔性多体动力学的共旋坐标法[J]. 力学季刊, 2017(2): 23 − 40.

    Shi Jiabei, Liu Zhuyong, Hong Jiazhen. The co-rotational formulation for flexible multibody dynamics [J]. Chinese Quarterly of Mechanics, 2017(2): 23 − 40. (in Chinese)
    [10] 王刚, 齐朝晖, 汪菁. 含铰接杆系结构几何非线性分析子结构方法[J]. 力学学报, 2014, 46(2): 273 − 283. doi:  10.6052/0459-1879-13-345

    Wang Gang, Qi Zhaohui, Wang Jing. Substructure methods of geometric nonlinear analysis for member structures with hinged supports [J]. Chinese Journal of Theoretical and Applied Mechani, 2014, 46(2): 273 − 283. (in Chinese) doi:  10.6052/0459-1879-13-345
    [11] 蔡松柏, 沈蒲生. 大转动平面梁有限元分析的共旋坐标法[J]. 工程力学, 2006, 23(增刊 1): 69 − 72. doi:  10.3969/j.issn.1000-4750.2006.z1.014

    Cai Songbai, Shen Pusheng. Co-rotational procedure for finite element analysis of plane beam element of large rotational displacement [J]. Engineering Mechanics, 2006, 23(Suppl 1): 69 − 72. (in Chinese) doi:  10.3969/j.issn.1000-4750.2006.z1.014
    [12] Battini J M. A non-linear corotational 4-node plane element [J]. Mechanics Research Communications, 2008, 35(6): 408 − 413. doi:  10.1016/j.mechrescom.2008.03.002
    [13] Li Z X. A stabilized co-rotational curved quadrilateral composite shell element [J]. International Journal for Numerical Methods in Engineering, 2011, 86(8): 975 − 999. doi:  10.1002/nme.3084
    [14] 周凌远, 李彤梅. 实体单元几何非线性CR列式有限元方法 [C]// 第18届全国结构工程学术会议论文集, 广州, 2009, Ⅰ: 562 − 566.

    Zhou Lingyuan, Li Tongmei. A co-rotational formulation for geometrically nonlinear finite element analysis of solid [C]// Proceedings of the 18th National Conference on Structural Engineering, Guangzhou, 2009, I: 562 − 566. (in Chinese)
    [15] 杜轲, 滕楠, 孙景江, 等. 基于共旋坐标和力插值纤维单元的RC框架结构连续倒塌构造方法[J]. 工程力学, 2019, 36(3): 105 − 114.

    Du Ke, Teng Nan, Sun Jingjiang, et al. A progressive collapse analytical model of RC frame structures based on corotational formulation for forced-based fiber elements [J]. Engineering Mechanics, 2019, 36(3): 105 − 114. (in Chinese)
    [16] 邓继华, 周福霖, 谭平. 圆钢管混凝土拱空间极限荷载计算方法研究[J]. 建筑结构学报, 2014, 35(11): 28 − 35.

    Deng Jihua, Zhou Fulin, Tan Ping. A study on method for calculating spatial ultimate load of circular CFST arch [J]. Journal of Building structure, 2014, 35(11): 28 − 35. (in Chinese)
    [17] 陈政清. 梁杆结构几何非线性有限元的数值实现方法[J]. 工程力学, 2014, 31(6): 42 − 52.

    Chen Zhengqing. Numerical implementation of geometrically nonlinear finite element method for beam structures [J]. Engineering Mechanics, 2014, 31(6): 42 − 52. (in Chinese)
    [18] Saafan S A. Theorefical analysis of suspension bridges [J]. Journal of the Structural Division, ASCE, 1966, 92(ST4): 1 − 11.
    [19] Nazmy A S. Three-dimensional nonlinear static analysis of cable-stayed bridges [J]. Computers & Structures, 1990, 34(2): 257 − 271.
    [20] 舒赣平, 刘伟, 陈绍礼. 半刚性钢框架的直接分析方法理论研究[J]. 建筑结构学报, 2014, 35(8): 142 − 150.

    Shu Ganping, Liu Wei, Chen Siulai. Theoretical research on direct analysis method for semi-rigid steel frames [J]. Journal of Building structure, 2014, 35(8): 142 − 150. (in Chinese)
    [21] 胡淑军, 王湛, 潘建荣. 基于截面组合法和截面弹簧刚度考虑跨内塑性铰的钢框架高等分析[J]. 工程力学, 2014, 31(2): 203 − 209.

    Hu Shujun, Wang Zhan, Pan Jianrong. Advanced analysis of steel frames using element with internal plastic hinge based on section assemblage concept and section spring stiffness [J]. Engineering Mechanics, 2014, 31(2): 203 − 209. (in Chinese)
    [22] 陈星烨, 邵旭东, 颜东煌. 大跨度斜拉桥地震动线性与非线性响应分析[J]. 湖南大学学报(自然科学版), 2002, 29(2): 106 − 111.

    Chen Xingye, Shao Xudong, Yan Donghuang. Long span bridge linear and nonlinear respond analysis while earthquake [J]. Journal of Hunan University (Natural Sciences Edition), 2002, 29(2): 106 − 111. (in Chinese)
    [23] 吕西林, 金国芳, 吴晓涵. 钢筋混凝土结构非线性有 限元理论与应用[M]. 上海: 同济大学出版社, 1997.

    Lü Xilin, Jin Guofang, Wu Xiaohan. Nonlinear finite element theory and application of reinforced concrete structure[M]. Shanghai: Tongji University Press, 1997. (in Chinese)
    [24] 蔡松柏, 沈蒲生. 关于非线性方程组求解技术[J]. 湖南大学学报(自然科学版), 2000, 27(3): 86 − 91.

    Cai Songbai, Shen Pusheng. On the methods solving nonlinear equations [J]. Journal of Hunan University (Natural Sciences Edition), 2000, 27(3): 86 − 91. (in Chinese)
    [25] Crisfield M A. A fast incremental/iterative solution procedure that handles“snap-through” [J]. Computers & Structures, 1981, 13(1/2/3): 55 − 62.
    [26] Lee S, Manuel F S, Rossow E C. Large deflections and stability of elastic frames [J]. Journal of Engineering Mechanics, 1968, 94(4): 521 − 547.
    [27] Remo Magalhaes de Souza. Force-based finite element for large displacement inelastic analysis of frames [D]. UC Berkeley: Department of Civil and Environmental Engineering, 2000.
    [28] 邓继华, 邵旭东. 基于场一致性的可带铰空间梁元几何非线性分析[J]. 工程力学, 2013, 30(1): 91 − 98.

    Deng Jihua, Shao Xudong. Geometric nonlinear analysis for 3D beam with hinge based on field consistency [J]. Engineering Mechanics, 2013, 30(1): 91 − 98. (in Chinese)
    [29] 陈至达. 杆、板、壳大变形理论 [M]. 北京: 科学出版社, 1996.

    Chen Zhida. Large deflection theory of truss, plate and shell [M]. Beijing: Science Press, 1996. (in Chinese)
  • [1] 陈朝晖, 杨帅, 杨永斌.  弹性膜结构几何非线性分析的刚体准则法 . 工程力学, 2020, 37(6): 246-256. doi: 10.6052/j.issn.1000-4750.2019.08.0488
    [2] 王刚, 齐朝晖, 孔宪超.  含机构位移起重机主副臂组合臂架结构几何非线性分析 . 工程力学, 2015, 32(7): 210-218. doi: 10.6052/j.issn.1000-4750.2013.12.1176
    [3] 陈政清.  梁杆结构几何非线性有限元的数值实现方法 . 工程力学, 2014, 31(6): 42-52. doi: 10.6052/j.issn.1000-4750.2013.05.ST08
    [4] 邓继华, 邵旭东.  钢管混凝土拱空间高效双非线性有限元分析 . 工程力学, 2014, 31(12): 89-95,111. doi: 10.6052/j.issn.1000-4750.2013.06.0565
    [5] 王振, 孙秦.  基于共旋三角形厚薄通用壳元的几何非线性分析 . 工程力学, 2014, 31(5): 27-33. doi: 10.6052/j.issn.1000-4750.2012.12.0958
    [6] 王庆, 姚竞争.  板壳结构非线性随机有限元法研究 . 工程力学, 2013, 30(12): 286-292. doi: 10.6052/j.issn.1000-4750.2012.06.0444
    [7] 苏静波, 范晓晨, 邵国建.  几何非线性扩展有限元法及其断裂力学应用 . 工程力学, 2013, 30(4): 42-46. doi: 10.6052/j.issn.1000-4750.2011.10.0702
    [8] 邓继华, 邵旭东.  基于混合单元法的预应力混凝土梁非线性分析 . 工程力学, 2013, 30(10): 171-177. doi: 10.6052/j.issn.1000-4750.2012.07.0494
    [9] 邓继华, 邵旭东.  基于场一致性的可带铰空间梁元几何非线性分析 . 工程力学, 2013, 30(1): 91-98. doi: 10.6052/j.issn.1000-4750.2011.07.0443
    [10] 邓继华, 邵旭东.  基于共旋坐标法的带刚臂平面梁元非线性分析 . 工程力学, 2012, 29(11): 143-151. doi: 10.6052/j.issn.1000-4750.2011.03.0123
    [11] 秦 剑, 黄克服, 张清东.  几何非线性样条有限元法 . 工程力学, 2011, 28(增刊I): 1-004.
    [12] 邓继华, 邵旭东, 邓潇潇.  四边形八节点共旋法平面单元的几何非线性分析 . 工程力学, 2011, 28(7): 6-012.
    [13] 蔡松柏, 沈蒲生, 胡柏学, 邓继华.  基于场一致性的2D四边形单元的共旋坐标法 . 工程力学, 2009, 26(12): 31-034.
    [14] 张年文, 童根树.  平面框架几何非线性分析的修正拉格朗日-协同转动联合法 . 工程力学, 2009, 26(8): 100-106,.
    [15] 蔡松柏, 沈蒲生.  大转动平面梁有限元分析的共旋坐标法 . 工程力学, 2006, 23(S1): 69-72,6.
    [16] 李世荣, 宋曦, 周又和.  弹性曲梁几何非线性精确模型及其数值解 . 工程力学, 2004, 21(2): 129-133.
    [17] 杨加明, 孙良新, 刘志和, 鲁宇明.  高阶剪切变形理论下两邻边铰支两邻边夹紧复合材料层板的几何非线性分析 . 工程力学, 2003, 20(4): 92-98.
    [18] 张少焱, 陈万吉.  基于几何非线性不协调元变分原理的圆柱壳非线性初始稳定性分析 . 工程力学, 2002, 19(1): 14-18.
    [19] 陈立群, 程昌钧, 张能辉.  具有几何和物理非线性粘弹性梁的混沌运动 . 工程力学, 2001, 18(1): 1-6.
    [20] 唐建民, 卓家寿, 何署廷.  索穹顶结构非线性分析的杆索有限元法 . 工程力学, 1998, 15(4): 34-42.
  • 加载中
图(7) / 表 (3)
计量
  • 文章访问数:  51
  • HTML全文浏览量:  21
  • PDF下载量:  16
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-01-07
  • 修回日期:  2020-07-16
  • 网络出版日期:  2020-11-06
  • 刊出日期:  2020-11-25

基于共旋法与稳定函数的几何非线性平面梁单元

doi: 10.6052/j.issn.1000-4750.2020.01.0012
    基金项目:  广东省自然科学基金项目(2015A030310141);中国博士后基金项目(2014M562154);国家自然科学基金项目(51478049);湖南省科技重大专项(2015GK1001-1);长沙理工大学土木工程优势特色重点学科创新性项目(2015ZDXK03);长沙理工大学研究生“实践创新与创业能力提升计划”项目(SJCX202001);长沙理工大学青年教师成长计划项目(2019QJCZ059)
    作者简介:

    谭建平(1996−),男,湖南人,硕士生,主要从事大跨径桥梁施工监控研究(E-mail: 834115536@qq.com)

    谭 平(1973−),男,湖南人,研究员,博士,博导,主要从事结构抗震和防灾减灾研究(E-mail: ptan@gzu.edu.cn)

    田仲初(1963−),男,湖南人,教授,博士,博导,主要从事大跨径桥梁施工监控研究(E-mail: 382525361@qq.com)

    通讯作者: 邓继华(1975−),男,湖南冷水江人,副教授,博士,硕导,主要从事桥梁与结构非线性研究(E-mail: jihuadeng@sina.com)

摘要: 建立一个准确、高效的几何非线性梁单元对于描述杆系结构的非线性行为至关重要。该文基于共旋坐标法和稳定函数提出了一种几何非线性平面梁单元。该单元在形成中把变形和刚体位移分开,局部坐标系内采用稳定函数以考虑单元P-δ效应的影响,从局部坐标系到结构坐标系的转换则采用共旋坐标法以及微分以考虑几何非线性,给出了几何非线性平面梁单元在结构坐标系下的全量平衡方程和切线刚度矩阵;在此基础上根据带铰梁端弯矩为零的受力特征,导出了能考虑梁端带铰的单元切线刚度矩阵表达式。通过多个典型算例验证了算法与程序的正确性、计算精度和效率。

English Abstract

邓继华, 谭建平, 谭平, 田仲初. 基于共旋法与稳定函数的几何非线性平面梁单元[J]. 工程力学, 2020, 37(11): 28-35. doi: 10.6052/j.issn.1000-4750.2020.01.0012
引用本文: 邓继华, 谭建平, 谭平, 田仲初. 基于共旋法与稳定函数的几何非线性平面梁单元[J]. 工程力学, 2020, 37(11): 28-35. doi: 10.6052/j.issn.1000-4750.2020.01.0012
Ji-hua DENG, Jian-ping TAN, Ping TAN, Zhong-chu TIAN. A GEOMETRIC NONLINEAR PLANE BEAM ELEMENT BASED ON COROTATIONAL FORMULATION AND ON STABILITY FUNCTIONS[J]. Engineering Mechanics, 2020, 37(11): 28-35. doi: 10.6052/j.issn.1000-4750.2020.01.0012
Citation: Ji-hua DENG, Jian-ping TAN, Ping TAN, Zhong-chu TIAN. A GEOMETRIC NONLINEAR PLANE BEAM ELEMENT BASED ON COROTATIONAL FORMULATION AND ON STABILITY FUNCTIONS[J]. Engineering Mechanics, 2020, 37(11): 28-35. doi: 10.6052/j.issn.1000-4750.2020.01.0012
  • 随着计算机软硬件技术的发展,基于非线性有限元开展结构破坏乃至失效倒塌分析已经成为可能[1-2],这其中研究开发出各种高效精确的非线性单元模型和算法是基础和前提[3-5]。对于几何非线性梁单元,许多学者从考虑应变的高阶项推导精确的单元切线刚度矩阵、杆端力与变形后单元坐标转换矩阵的精确计算等入手,进行了卓有成效的研究[4, 6]。在各种几何非线性有限元列式中,共旋坐标(CR)法认为结构的几何非线性主要由结构大位移引起,因此在分析中每一个单元都建立一个共旋坐标系(局部坐标系),该坐标系随单元位移发生平移和转动,该坐标系下单元的纯变形可以通过扣除单元刚体平移和转动而得到,相应可求得单元切线刚度矩阵和抗力,再通过变换的坐标转换矩阵将其转换到结构坐标系下,此过程中几何非线性效应被自动计入[7-9]。研究表明,基于共旋坐标法研究几何非线性,不仅具有列式简单、计算准确,以及适用于杆、梁、平面、壳、实体单元等各类型单元的优点[10-14],而且能将几何与材料非线性脱藕,即材料非线性在局部坐标系内考虑,几何非线性则通过局部坐标系与结构坐标系之间的转换予以实现[15-16]。对于梁单元而言,除了大位移能引起几何非线性以外,梁-柱效应也是引起几何非线性的一个重要因素[17],这其中采用能考虑构件轴力与弯矩相互作用的稳定函数[18-19]来考虑梁柱效应是一种最精确的方法[20-21],实际工程分析中常采用的几何刚度矩阵实质是将稳定函数法表达的单元刚度矩阵用级数展开,取其一阶近似而得到[22]。但应指出的是,以往基于稳定函数法的研究采用的有限元列式一般为T.L或U.L列式,笔者尚未见到将稳定函数法与共旋坐标(CR)法结合起来进行几何非线性平面梁单元的研究。因此,本文基于共旋坐标法,在局部坐标系下分离单元刚体位移和变形,采用稳定函数考虑梁柱效应;再基于结构坐标系与局部坐标系下节点力之间及节点位移之间的总量关系及微分导出的增量关系,获得平面梁单元在结构坐标系下的全量平衡方程和切线刚度矩阵;在此基础上根据带铰梁端弯矩为零的受力特征,导出了能考虑梁端带铰的单元切线刚度矩阵。最后通过数个经典算例对本文算法及程序进行了较全面和严格的检验。

    • 图1(a)、图1(b)所示分别为初始时刻和计算$t$时刻的平面梁单元,$XY$为结构坐标系,$xy$为单元的共旋坐标系,它具有随单元变形而转动的特点,在运动中始终以节点$i$为原点,以节点$i$$j$的连线方向为$x$轴,$y$轴则由$x$轴按逆时钟方向旋转90°形成。

      图  1  变形前后平面梁单元

      Figure 1.  Plane beam element before and after deformation

      初始时刻设单元在结构坐标系下水平及竖向的投影长度为${x^0}$${y^0}$,在计算$t$时刻结构坐标系中的位移向量${{d}}={[\begin{array}{*{20}{c}} {{u_i}}&{{v_i}}&{{\theta _i}}&{{u_j}}&{{v_j}}&{{\theta _j}} \end{array}]^{\rm{T}}}$,在共旋坐标系中的位移向量${{{d}}^{\rm{L}}}=[ \begin{array}{*{20}{c}} {u_i^{\rm{L}}}&{v_i^{\rm{L}}}&{\theta _i^{\rm{L}}}&{u_j^{\rm{L}}} &{v_j^{\rm{L}}}\end{array} $ $ {\begin{array}{*{20}{c}} {\theta _j^{\rm{L}}} \end{array}} ]^{\rm{T}}$,联立图1(a)、图1(b),有:

      $$ \begin{split} & {x^t}={x^0} + {u_j} - {u_i},\\& {y^t}={y^0} + {v_j} - {v_i}, \\& {\alpha ^t}=\arctan \left(\frac{{{y^t}}}{{{x^t}}}\right) \end{split} $$ (1)

      平面梁单元在共旋坐标系中的节点位移为:

      $$ \begin{split} & u_i^{\rm{L}}=v_i^{\rm{L}}=v_j^{\rm{L}}=0,\;\;\;\;\;\;\;\;\; u_j^{\rm{L}}={l^t} - {l^0}, \\& \theta _i^{\rm{L}}={\theta _i} - ({\alpha ^t} - {\alpha ^0}),\;\;\;\;\;\; \theta _j^{\rm{L}}={\theta _j} - ({\alpha ^t} - {\alpha ^0}) \end{split} $$ (2)

      式中,${l^0}$${l^t}$分别为变形前后梁单元长度,且有:

      $$ \begin{split} & {l^0}=\sqrt {{{( {{x^0}} )}^2} + {{( {{y^0}} )}^2}}, \\& {l^t}=\sqrt {{{( {{x^0} + {u_j} - {u_i}} )}^2} + {{( {{y^0} + {v_j} - {v_i}} )}^2}} \end{split} $$ (3)

      对式(2)微分且联立式(1)和式(3),可将共旋坐标系下位移微分${\rm{\delta }}{{{d}}^{\rm{L}}}=[ {\begin{array}{*{20}{c}} {{\rm{\delta }}u_i^{\rm{L}}}&{{\rm{\delta }}v_i^{\rm{L}}}&{{\rm{\delta }}\theta _i^{\rm{L}}} &{{\rm{\delta }}u_j^{\rm{L}}}&{{\rm{\delta }}v_j^{\rm{L}}}\end{array}} $ ${ {\begin{array}{*{20}{c}} {{\rm{\delta }}\theta _j^{\rm{L}}} \end{array}} ]^{\rm{T}}}$用结构坐标系下位移微分${\rm{\delta }}{{d}}=[ {\begin{array}{*{20}{c}} {{\rm{\delta }}{u_i}}&{{\rm{\delta }}{v_i}} \end{array}}$${ {\begin{array}{*{20}{c}} {{\rm{\delta }}{\theta _i}}&{{\rm{\delta }}{u_j}}&{{\rm{\delta }}{v_j}}&{{\rm{\delta }}{\theta _j}} \end{array}} ]^{\rm{T}}} $表示为:

      $${\rm{\delta }}{{{d}}^{\rm{L}}}={{T}} \cdot {\rm{\delta }}{{d}}$$ (4)

      其中:

      $${{T}}=\left[ {\begin{array}{*{20}{c}} 0&0&0&0&0&0 \\ 0&0&0&0&0&0 \\ { - \dfrac{{si}}{{{l^t}}}}&{\dfrac{{co}}{{{l^t}}}}&1&{\dfrac{{si}}{{{l^t}}}}&{ - \dfrac{{co}}{{{l^t}}}}&0 \\ { - co}&{ - si}&0&{co}&{si}&0 \\ 0&0&0&0&0&0 \\ { - \dfrac{{si}}{{{l^t}}}}&{\dfrac{{co}}{{{l^t}}}}&0&{\dfrac{{si}}{{{l^t}}}}&{ - \dfrac{{co}}{{{l^t}}}}&1 \end{array}} \right]$$

      式中:$si=\sin {\alpha ^t},\;{co=\cos \alpha } ^t$

      图2(a)、图2(b)所示分别为梁单元及微段的受力平衡情况。

      图  2  梁单元及微段

      Figure 2.  beam element and micro-segment

      当单元轴力为拉力(以拉为正)时,由图2(b)知:

      $${\rm{d}}M= - Q{\rm{d}}x - N{\rm{d}}y\qquad\qquad\;\; $$ (5)

      微分式(5),有:

      $$\frac{{{{\rm{d}}^2}M}}{{{\rm{d}}{x^2}}}= - N\frac{{{{\rm{d}}^2}y}}{{{\rm{d}}{x^2}}}\qquad\qquad\;\;\;\;\;\;\; $$ (6)

      在初等梁理论中因$M= - EI\dfrac{{{{\rm{d}}^2}y}}{{{\rm{d}}{x^2}}}$,故有:

      $$\frac{{{{\rm{d}}^4}y}}{{{\rm{d}}{x^4}}} - \frac{N}{{EI}}\frac{{{{\rm{d}}^2}y}}{{{\rm{d}}{x^2}}}=0\qquad\qquad\;\;\; $$ (7)

      要求解式(7)所示的常系数四阶微分方程,须利用以下位移和杆端力的边界条件,即:

      $$ \begin{split} & y\left| {_{x=0}} \right.=0,\;\;\;\;y\left| {_{x={l^t}}} \right.=0 ,\\& {y'}\left| {_{x=0}} \right.=\theta _i^{\rm{L}},\;\;\;\;y'\left| {_{x={l^t}}} \right.=\theta _j^{\rm{L}} ,\\& M\left| {_{x=0}} \right.={M_i},\;\;\;\;M\left| {_{x={l^t}}} \right.= - {M_j} \end{split} $$ (8)

      利用式(8)所示的边界条件解得挠度$y$后,再进行相关运算可得到杆端弯矩${M_i}$${M_j}$的计算式为:

      $$ \begin{split} & {M_i}=sk\theta _i^{\rm{L}} + sck\theta _j^{\rm{L}}, \\& {M_j}=sck\theta _i^{\rm{L}} + sk\theta _j^{\rm{L}} \end{split} \qquad\quad\;\;$$ (9)

      式中各参数$s$$c$$k$$\omega $的计算式为[18]

      $$ \begin{split} & s=\frac{{\omega (2\omega {\rm{ch}} 2\omega - {\rm{sh}} 2\omega )}}{{1 - {\rm{ch}} 2\omega + \omega {\rm{sh}} 2\omega }},\\& c=\frac{{{\rm{sh}} 2\omega - 2\omega }}{{2\omega {\rm{ch}} 2\omega - {\rm{sh}} 2\omega }} ,\\& k=\frac{{EI}}{{{l^t}}},\\& \omega =\frac{{{l^t}}}{2}\sqrt {\frac{N}{{EI}}} \end{split}\quad\;\; $$ (10)

      当单元轴力为压力时,式(5)就变为下式:

      $${\rm{d}}M= - Q{\rm{d}}x + N{\rm{d}}y\qquad\;\;\;\;$$ (11)

      对式(11)进行与单元受拉时同样的推导,可得到与式(9)相同的表达式,只是各参数$s$$c$$k$$\omega $的计算式变成以下形式:

      $$ \begin{split} & s=\frac{{\omega (1 - 2\omega {\rm{ctg}} 2\omega )}}{{{\rm{tg}} \omega - \omega }},\\& c=\frac{{\sin 2\omega - 2\omega }}{{2\omega \cos 2\omega - \sin 2\omega }} ,\\& k=\frac{{EI}}{{{l^t}}},\\& \omega =\frac{{{l^t}}}{2}\sqrt {\frac{{ - N}}{{EI}}} \end{split} $$ (12)

      应当指出的是,当单元轴力为0时,由洛比塔法则可得到$s$=4和$c$=0.5。

      设单元在共旋坐标系下的杆端力向量为${{f}}=$${[ {\begin{array}{*{20}{c}} {{N_i}}&{{Q_i}}&{{M_i}}&{{N_j}}&{{Q_j}}&{{M_j}} \end{array}} ]^{\rm{T}}}$,依对应及平衡关系显然有:

      $$ \begin{split} & {M_i}={M_i}, \\& {M_j}={M_j} ,\\& {Q_i}= - \frac{{({M_1} + {M_2})}}{{{l^t}}}= - s(1 + c)k(\theta _i^{\rm{L}} + \theta _j^{\rm{L}}) ,\\& {Q_j}=s(1 + c)k(\theta _i^{\rm{L}} + \theta _j^{\rm{L}}) ,\\& {N_i}=N=\frac{{EA}}{{{l^0}}}u_j^{\rm{L}} ,\\& {N_j}= - N= - \frac{{EA}}{{{l^0}}}u_j^{\rm{L}} \end{split} $$ (13)

      为方便计,可将式(13)所示的共旋坐标系下杆端力${{f}}$与杆端位移${{{d}}^{\rm{L}}}$之间的关系用矩阵形式表示为:

      $${{f}}={{{k}}_{\rm{L}}} \cdot {{{d}}^{\rm{L}}}$$ (14)

      式中: $ {{{k}}_{\rm{L}}}=$

      $$ \left[ \!\!{\begin{array}{*{20}{c}} {0}&0&0&{\dfrac{{EA}}{{{l^0}}}}&0&0\\ 0&0&{ - s\left( {1 + c} \right)k}&0&0&{ - s\left( {1 + c} \right)k}\\ 0&0&{sk}&0&0&{sck}\\ 0&0&0&{ - \dfrac{{EA}}{{{l^0}}}}&0&0\\ 0&0&{s\left( {1 + c} \right)k}&0&0&{s\left( {1 + c} \right)k}\\ 0&0&{sck}&0&0&{sk} \end{array}} \!\!\right] $$ (15)

      设平面梁单元在结构坐标系的节点力向量为${{F}}={[ {\begin{array}{*{20}{c}} {F_i^X}&{F_i^Y}&{M_i^Z}&{F_j^X}&{F_j^Y}&M _j^Z \end{array}}]^{\rm{T}}}$,由场一致性原则有[11, 16]

      $${{F}}={{{t}}^{\rm{T}}} \cdot {{f}}$$ (16)

      式中:

      ${{t}}=\left[ {\begin{array}{*{20}{c}} {co}&{si}&0&0&0&0 \\ { - si}&{co}&0&0&0&0 \\ 0&0&1&0&0&0 \\ 0&0&0&{co}&{si}&0 \\ 0&0&0&{ - si}&{co}&0 \\ 0&0&0&0&0&1 \end{array}} \right]$

      微分式(16)可得:

      $${\rm{\delta }}{{F}}={\rm{\delta }}{{{t}}^{\rm{T}}} \cdot {{f}} + {{{t}}^{\rm{T}}} \cdot {\rm{\delta }}{{f}}$$ (17)

      为便于推导,可将${\rm{\delta }}{{{t}}^{\rm{T}}} \cdot {{f}}$表示成:

      $${\rm{\delta }}{{{t}}^{\rm{T}}} \cdot {{f}}=\left\{ {\begin{array}{*{20}{c}} {{\rm{\delta }}{{co}}{N_i} - {\rm{\delta }}{{si}}{Q_i}} \\ {{\rm{\delta }}{{si}}{N_i} + {\rm{\delta }}{{co}}{Q_i}} \\ 0 \\ {{\rm{\delta }}{{co}}{N_j} - {\rm{\delta }}{{si}}{Q_j}} \\ {{\rm{\delta }}{{si}}{N_j} + {\rm{\delta }}{{co}}{Q_j}} \\ 0 \end{array}} \right\}$$ (18)

      其中:

      $$\begin{array}{l} {\rm{\delta }}{{co}}= - \dfrac{{si}}{{{}^tl}} \cdot {{AA}} \cdot {\rm{\delta }}{{d}} ,\\ {\rm{\delta }}{{si}}=\dfrac{{co}}{{{}^tl}} \cdot {{AA}} \cdot {\rm{\delta }}{{d}}, \\ {{AA}}=[ {\begin{array}{*{20}{c}} {si}&{ - co}&0&{ - si}&{co}&0 \end{array}} ] \end{array} $$

      为便于表述,将式(18)改写为:

      $${\rm{\delta }}{{{t}}^{\rm{T}}} \cdot {{f}}={{{K}}_{\rm{n}}} \cdot {\rm{\delta }}{{d}}$$ (19)

      式中:

      $${{{K}}_{\rm{n}}}=\left[ {\begin{array}{*{20}{c}} { - \left( {\dfrac{{si}}{{{}^tl}}{N_i} + \dfrac{{co}}{{{}^tl}}{Q_i}} \right) \cdot {{AA}}} \\ {\left( {\dfrac{{co}}{{{}^tl}}{N_i} - \dfrac{{si}}{{{}^tl}}{Q_i}} \right) \cdot {{AA}}} \\ 0 \\ { - \left( {\dfrac{{si}}{{{}^tl}}{N_j} + \dfrac{{co}}{{{}^tl}}{Q_j}} \right) \cdot {{AA}}} \\ {\left( {\dfrac{{co}}{{{}^tl}}{N_j} - \dfrac{{si}}{{{}^tl}}{Q_j}} \right) \cdot {{AA}}} \\ 0 \end{array}} \right]$$

      式中,${N_i}$${Q_i}$${N_j}$${Q_j}$的值通过式(13)计算。

      为将${{{t}}^{\rm{T}}} \cdot {\rm{\delta }}{{f}}$$ \delta {{d}} $表示,可按以下步骤进行:微分式(3)可得到$ {\rm{\delta}} {{{l}}^{{t}}}$${\rm{\delta }}{{d}}$表达的式子,由式(2)知${\rm{\delta }}{{{l}}^{{t}}}{\rm{=\delta }}{{u}}_j^{\rm{L}}$,至此很容易由式(13)得到${\rm{\delta }}{{N}}$${\rm{\delta }}{{d}}$表达的计算式,在此基础上,根据单元受拉还是受压分别微分式(10)或式(12),可得到${\rm{\delta }}{{s}}$${\rm{\delta }}{{c}}$${\rm{\delta }}{{k}}$${\rm{\delta }}{{\omega }}$分别用${\rm{\delta }}{{d}}$表达的计算式,再微分式(14)并联立${\rm{\delta }}{{{l}}^{{t}}}{\rm{=\delta }}{{u}}_j^{\rm{L}}$及式(4),可得到${\rm{\delta }}{{f}}$${\rm{\delta }}{{d}}$表达的计算式为:

      $${\rm{\delta }}{{f}}={{{k}}^{\rm{T}}}{\rm{\delta }}{{d}}$$ (20)

      ${{{k}}^{\rm{T}}}$的具体表达式较繁琐,限于篇幅此处未列出。

      联立式(19)和式(20),可将式(17)写成:

      $$ {\rm{\delta }}{{F}}=( {{{{K}}_{\rm{n}}}{\rm{ + }}{{{t}}^{\rm{T}}}{{{k}}^{\rm{T}}}} ){\rm{\delta }}{{d}} $$ (21)

      因此,本文导出的结构坐标系下平面梁单元考虑几何非线性的切线刚度矩阵${{K}}$为:

      $${{K}}={{{K}}_{\rm{n}}}{\rm{ + }}{{{t}}^{\rm{T}}}{{{k}}^{\rm{T}}}$$ (22)
    • 假如铰在梁元的$i$端,此时有${M_i}=0$,由式(9)知,由于$s$$k$并不为0,故有:

      $$\theta _i^{\rm{L}}= - c\theta _j^{\rm{L}}$$ (23)

      由式(23)可知在梁元$i$端有铰的情况下,$\theta _i^{\rm{L}}$并不是一个独立的位移变量,联立式(23)和式(14)可消去$\theta _i^{\rm{L}}$,这种情况下${{{k}}_{\rm{L}}}$的表达式变成:

      $$ \begin{split} & {{{k}}_{\rm{L}}}=\\&\left[ {\begin{array}{*{20}{c}} {0}&0&0&{\dfrac{{EA}}{{{l^0}}}}&0&0\\ 0&0&{ - s( {1 + c} )k}&0&0&{ - s( {1 + c} )k}\\ 0&0&0&0&0&0\\ 0&0&0&{ - \dfrac{{EA}}{{{l^0}}}}&0&0\\ 0&0&{s( {1 + c} )k}&0&0&{s( {1 + c} )k}\\ 0&0&0&0&0&{sk( {1 - {c^2}} )} \end{array}} \right] \end{split} $$ (24)

      进行同样的推导,不难得到$i$端带铰时式(22)中的${{K}}$的具体表达式,当铰在梁元的$j$端时,同样可推导得出。

    • 非线性方程组的增量法求解主要有荷载增量法和位移增量法[23]。理论上而言,荷载增量法只能计算出荷载-位移曲线的上升段,而位移增量法能计算荷载-位移曲线的下降段,从而能获得极限荷载值与位移值;但在实际计算中发现,由于在荷载-位移曲线的顶点附近结构切线刚度接近奇异,即使采用位移增量法,非线性计算也很难收敛。鉴于此,很多文献都大力推荐弧长增量法,文献[24]对传统弧长增量法[25]研究后认为,该方法能顺利进行的前提是必须对结构的特性有深刻的了解并进行复杂的试算,显然这会给分析工作带来极大困难和不便;因此,该文献在传统弧长增量法的基础上提出一种如图3所示的基于牛顿-拉菲逊法的投影增量法,该方法能有效克服传统弧长增量法的上述缺点,且该方法收敛速度要稍快,计算量也略小,本文采用该方法进行非线性方程组的求解。

      本文非线性计算流程如下:

      1)根据式(2)由上一荷载步末单元$i$$j$节点在结构坐标系下的总位移向量${{d}}$求出共旋坐标系下的节点位移向量$ {{{d}}^{\rm{L}}}$

      2)由式(13)中最后一行单元轴力N的计算式得到N,再根据轴力N是受拉还是受压分别基于式(10)或式(12)计算参数$s$$c$$k$$\omega $,进而根据式(9)得到单元弯矩${M_i}$${M_j}$

      3)根据式(13)可得到单元在共旋坐标系下的杆端力向量为f

      4)由式(4)和式(16)分别计算矩阵Tt

      5)由式(19)计算Kn,再由式(22)得到单元在结构坐标系下的切线刚度矩阵K

      6)由式(16)得到单元在结构坐标系下的节点抗力向量F

      7)对所有单元重复步骤1)~步骤6),基于常规的“对号入座”原则叠加形成结构总刚矩阵和总抗力矩阵;

      8)由总抗力矩阵与总荷载矩阵之差形成不平衡力矩阵作为非线性平衡方程的右端项,解方程得到增量位移矩阵$ \Delta {{d}}$,与前述总位移向量d叠加形成新的总位移向量d

      9)判断是否收敛,如是,转至下一个荷载步;如否,进行下一轮迭代计算。

      图  3  投影增量法示意

      Figure 3.  Projection increment method

    • 采用3个经典算例对本文算法及程序从正确性、计算效率及精度方面进行较全面和严格的检验。

      算例1. 如图4所示为Lee’s框架[26],材料及几何参数如下:线弹性模量E=7060.8 kN/m2,梁和柱具有同样的截面,且截面面积为6 cm2,惯性矩为2 cm4,对结构进行几何非线性分析,分析中柱均分为8个单元,荷载作用点左边梁均分为2个单元,右边梁均分为6个单元,图5(a)所示为框架的初始构型以及不同变形时刻的构型,图5(b)图5(c)所示分别为荷载作用点处荷载-水平位移及荷载-竖向位移曲线,可看出与相关文献[27]的计算结果很吻合(该文献未提供荷载-水平位移曲线)。

      图  4  Lee’s框架及荷载 /cm

      Figure 4.  Lee’s frame geometry and loading

      图  5  Lee’s框架

      Figure 5.  Lee’s frame

      算例2 . 如图6所示悬臂梁,集中荷载$P$=35 kN作用在自由端,梁跨为$L$=100 m、弯曲刚度为$EI$=3.5×104 kN·m2、梁端水平位移为$u$、垂直位移为$w$、转角为$\theta $,对本算例按以下两种方法求解:

      图  6  端部承受集中荷载的悬臂梁

      Figure 6.  Cantilever subjected to concentrated load at free end

      方法1. 采用本文算法研制的程序进行计算;

      方法2. 按文献[11]和文献[28](文献[28]就是将文献[11]的算法从平面梁延伸到空间梁)的方法进行计算,即共旋坐标系下的杆端力就由节点位移(去除了刚体位移)与小位移线弹性刚度矩阵相乘得到。

      为使计算结果具有可比性,两个程序采用的方程组求解方法、数据精度及收敛精度标准等都完全相同,方法2中杆端力也是基于全量计算,详见文献[11]及文献[28]。

      两种方法得到的计算结果及比较见表1表2,解析解见文献[29]。可以看出,无论是从非线性计算能力、精度还是稳定性方面,方法1均优于方法2(方法2在划分为1个单元及4个荷载步情况下还出现错误的收敛解);还可以看出,由于方法1与方法2中杆端力的计算均是基于全量计算,所以在单元数划分固定的前提下,计算精度与荷载步数无关。

      表 1  悬臂梁自由端的位移

      Table 1.  Displacement of cantilever at free end

      荷载步数方法1个单元2个单元
      $u$$w$$\theta $$u$$w$$\theta $
      3152.33587.9181.45053.89383.4981.435
      256.76184.7781.478
      4152.33587.9181.45053.89383.4981.435
      2136.0593.2821.803
      5152.33587.9181.45053.89383.4981.435
      256.76184.7781.478
      6152.33587.9181.45053.89383.4981.435
      263.94493.2821.803
      7152.33587.9181.45053.89383.4981.435
      263.94493.2821.803
      解析解:$u$=55.5 $w$=81.06 $\theta $=1.430
      注:短横线表示该单元划分及荷载步条件下非线性计算失败。

      表 2  计算误差

      Table 2.  Calculation error /(%)

      方法1个单元2个单元
      $u$$w$$\theta $$u$$w$$\theta $
      15.708.461.382.903.010.35
      215.2115.0826.052.274.593.31

      算例3. 如图7(a)所示,一对拉力2$P$作用于铰接方棱形框架的对角点,各项几何、材料参数及位移变量详见图7,EI为杆件的抗弯刚度。

      图  7  对角点受拉力作用时的铰接方棱形框架

      Figure 7.  Diamond-shaped frame under a pare of opposite concentrated tensions in opposite angles

      为验证本文所导出的带铰梁单元模型,基于荷载及结构对称的特征,取图7(b)所示的1/2结构进行计算,按每根杆件分成1、2、3和4个单元(中间铰左、右两个相邻单元类型选用本文推导的带铰单元),在荷载参数${{P{L^2}} / {EI}}$=10时所得的计算结果如表3所示,其中解析解见文献[29]。从表3可看出,计算结果收敛非常好,如按工程精度误差不超过5%的要求,则每根杆件划分成3个单元就已经足够。

      表 3  铰接方棱形框架位移

      Table 3.  Displacement of diamond-shaped frame

      单元数$u/L$$w/L$${\theta _0}$最大误差/(%)
      10.46180.26247.6
      20.46190.25071.39976.9
      30.46520.24761.46382.6
      40.46510.24591.48161.5
      解析解0.46600.24381.5035
      注:在每根杆件划分成1个单元时,${\theta _0}$=0显然失真。
    • 本文在局部坐标系下采用稳定函数考虑P-δ效应,再基于结构坐标系与局部坐标系下节点力之间及节点位移之间的总量关系及微分导出的增量关系,获得平面梁单元在结构坐标系下的切线刚度矩阵,同时还获得单元杆端力的全量算法;在此基础上,根据带铰梁端弯矩为零的受力特征,导出了能考虑梁端带铰的单元切线刚度矩阵表达式。多个算例表明本文算法是正确的,相对于已有文献算法,本文算法在非线性计算能力、精度以及稳定性方面均占有一定优势,可用于平面杆系结构的几何非线性分析。

参考文献 (29)

目录

    /

    返回文章
    返回