留言板

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

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

功能梯度梁静力与动力分析的格林函数法

邓先琪 苏成 马海涛

邓先琪, 苏成, 马海涛. 功能梯度梁静力与动力分析的格林函数法[J]. 工程力学. doi: 10.6052/j.issn.1000-4750.2019.10.0615
引用本文: 邓先琪, 苏成, 马海涛. 功能梯度梁静力与动力分析的格林函数法[J]. 工程力学. doi: 10.6052/j.issn.1000-4750.2019.10.0615
Xian-qi DENG, Cheng SU, Hai-tao MA. GREEN’S FUNCTION METHOD FOR STATIC AND DYNAMIC ANALYSIS OF FUNCTIONALLY GRADED BEAMS[J]. Engineering Mechanics. doi: 10.6052/j.issn.1000-4750.2019.10.0615
Citation: Xian-qi DENG, Cheng SU, Hai-tao MA. GREEN’S FUNCTION METHOD FOR STATIC AND DYNAMIC ANALYSIS OF FUNCTIONALLY GRADED BEAMS[J]. Engineering Mechanics. doi: 10.6052/j.issn.1000-4750.2019.10.0615

功能梯度梁静力与动力分析的格林函数法

doi: 10.6052/j.issn.1000-4750.2019.10.0615
基金项目: 国家自然科学基金项目(51678252);广州市科学研究计划重点项目(201804020069)
详细信息
    作者简介:

    邓先琪(1989−),男,广东人,博士生,主要从事结构不确定性问题研究(E-mail: dengxqi@126.com)

    马海涛(1962−),男,辽宁人,教授,博士,博导,主要从事计算力学、结构分析与优化理论与算法研究(E-mail: htma@gzhu.edu.cn)

    通讯作者: 苏 成(1968−),男,广东人,教授,博士,博导,主要从事结构随机振动和计算力学研究(E-mail: cvchsu@scut.edu.cn)
  • 中图分类号: O342

GREEN’S FUNCTION METHOD FOR STATIC AND DYNAMIC ANALYSIS OF FUNCTIONALLY GRADED BEAMS

  • 摘要: 功能梯度梁静动态响应的数值分析方法一般局限于有限元法,存在有限元法的固有缺点,有必要发展新的数值求解方法。将功能梯度梁静力分析的控制微分方程转化为与匀质材料梁静力分析控制微分方程相一致的形式,并利用匀质材料梁静力问题的格林函数,开展功能梯度梁的静力分析。在此基础上,进一步推导获得功能梯度梁的柔度矩阵,据此建立功能梯度梁的运动方程,开展功能梯度梁的动力特性分析和动力响应分析。数值算例表明,采用格林函数法可以高效准确地分析功能梯度梁的静力响应与动力行为,验证了方法的计算精度与效率。
  • 图  1  与功能梯度梁等效的匀质梁

    Figure  1.  A uniform beam equivalently to the functionally graded beam

    图  2  一端固支、一端简支功能梯度梁

    Figure  2.  A clamped-simply supported functionally graded beam

    图  3  功能梯度梁跨中点C的轴向与横向位移

    Figure  3.  Axial and transverse displacement of node C at midspan of the functionally graded beam

    图  4  功能梯度梁固定端点A的轴力和弯矩

    Figure  4.  Axial force and bending moment of node A at fixed end of the functionally graded beam

    图  5  功能梯度梁前三阶横向自由振动固有频率

    Figure  5.  The first three natural frequencies of functionally graded beam’s transverse free vibration

    图  6  荷载-时间曲线

    Figure  6.  Load-time curve

    图  7  功能梯度梁跨中点C横向位移、速度和加速度

    Figure  7.  Transverse displacement, velocity, acceleration of node C at midspan of the functionally graded beam

    表  1  功能梯度梁跨中点C的轴向与横向位移

    Table  1.   Axial and transverse displacement of node C at midspan of the functionally graded beam

    位移方法子段数n,离散单元数Ne
    4102050100
    轴向位移/
    (×10−5 m)
    解析解 1.328
    有限元法 1.314 1.326 1.327 1.328 1.328
    格林函数法
    (d1=d2=0.2 m)
    1.349 1.331 1.329 1.328 1.328
    格林函数法
    (d1=d2=0.4 m)
    1.349 1.331 1.329 1.328 1.328
    格林函数法
    (d1=0.2 m, d2=0.2 m)
    1.349 1.331 1.329 1.328 1.328
    横向位移/
    (×10−3 m)
    解析解 −2.447
    有限元法 −2.403 −2.439 −2.445 −2.447 −2.447
    格林函数法
    (d1=d2=0.2 m)
    −2.489 −2.451 −2.448 −2.447 −2.447
    格林函数法
    (d1=d2=0.4 m)
    −2.489 −2.451 −2.448 −2.447 −2.447
    格林函数法
    (d1=0.2 m, d2=0.2 m)
    −2.489 −2.451 −2.448 −2.447 −2.447
    下载: 导出CSV

    表  2  功能梯度梁固定端点A的轴力与弯矩

    Table  2.   Axial force and bending moment of node A at fixed end of the functionally graded beam

    内力方法子段数n,离散单元数Ne
    4102050100
    轴力/kN 解析解 3.166
    有限元法 3.181 3.168 3.166 3.166 3.166
    格林函数法
    (d1=d2=0.2 m)
    3.066 3.149 3.161 3.165 3.165
    格林函数法
    (d1=d2=0.4 m)
    3.066 3.149 3.161 3.165 3.165
    格林函数法
    (d1=0.2 m, d2=0.2 m)
    3.066 3.149 3.161 3.165 3.165
    弯矩/(kN$ \cdot $m) 解析解 −0.625
    有限元法 −0.634 −0.627 −0.626 −0.625 −0.625
    格林函数法
    (d1=d2=0.2 m)
    −0.677 −0.634 −0.628 −0.626 −0.625
    格林函数法
    (d1=d2=0.4 m)
    −0.677 −0.634 −0.628 −0.626 −0.625
    格林函数法
    (d1=0.2 m, d2=0.2 m)
    −0.677 −0.634 −0.628 −0.626 −0.625
    下载: 导出CSV

    表  3  功能梯度梁前三阶横向自由振动固有频率

    Table  3.   The first three natural frequencies of functionally graded beam’s transverse free vibration

    频率方法子段数n,离散单元数Ne
    1020304050200
    第一阶频率/
    (×103 rad/s)
    有限元法1.3201.3231.3241.3241.3241.324
    格林函数法
    (m=4)
    1.3221.3231.3231.3231.3231.323
    格林函数法
    (m=6)
    1.3231.3241.3241.3241.3241.324
    格林函数法
    (m=8)
    1.3231.3241.3241.3241.3241.324
    格林函数法
    (m=10)
    1.3221.3241.3241.3241.3241.324
    第二阶频率/
    (×103 rad/s)
    有限元法4.2754.3074.3124.3144.3154.317
    格林函数法
    (m=4)
    4.2984.3034.3034.3044.3044.304
    格林函数法
    (m=6)
    4.3084.3144.3154.3154.3164.316
    格林函数法
    (m=8)
    4.3084.3154.3164.3174.3174.317
    格林函数法
    (m=10)
    4.3094.3154.3164.3174.3174.317
    第三阶频率/
    (×103 rad/s)
    有限元法8.8538.9929.0189.0279.0319.038
    格林函数法
    (m=4)
    8.8858.8978.8988.8998.8998.899
    格林函数法
    (m=6)
    8.9898.9958.9948.9948.9948.993
    格林函数法
    (m=8)
    9.0119.0329.0349.0359.0359.036
    格林函数法
    (m=10)
    9.0109.0339.0369.0379.0379.038
    下载: 导出CSV
  • [1] Koizumi M. FGM activities in Japan [J]. Composites Part B: Engineering, 1997, 28(1/2): 1 − 4. doi:  10.1016/S1359-8368(96)00016-9
    [2] 毛丽娟, 马连生. 非均匀热载荷作用下功能梯度梁的非线性静态响应[J]. 工程力学, 2017, 34(6): 1 − 8.

    Mao Lijuan, Ma Liansheng. Nonlinear static responses of FGM beams under non-uniform thermal loading [J]. Engineering Mechanics, 2017, 34(6): 1 − 8. (in Chinese)
    [3] 何昊南, 于开平. 考虑热对材料参数影响的FGM梁热后屈曲特性研究[J]. 工程力学, 2019, 36(4): 52 − 61.

    He Haonan, Yu Kaiping. Thermal post-buckling analysis of FGM beams considering the heat effect on materials [J]. Engineering Mechanics, 2019, 36(4): 52 − 61. (in Chinese)
    [4] Celebi K, Yarimpabuc D, Tutuncu N. Free vibration analysis of functionally graded beams using complementary functions method [J]. Archive of Applied Mechanics, 2018, 88(5): 729 − 739. doi:  10.1007/s00419-017-1338-6
    [5] Cao D X, Gao Y H, Yao M H, et al. Free vibration of axially functionally graded beams using the asymptotic development method [J]. Engineering Structures, 2018, 173: 442 − 448. doi:  10.1016/j.engstruct.2018.06.111
    [6] Birman V, Byrd L W. Modeling and analysis of functionally graded materials and structures [J]. Applied Mechanics Reviews, 2007, 60(5): 195 − 216. doi:  10.1115/1.2777164
    [7] 仲政, 吴林志, 陈伟球. 功能梯度材料与结构的若干力学问题研究进展[J]. 力学进展, 2010, 40(5): 528 − 541.

    Zhong Zheng, Wu Linzhi, Cheng Weiqiu. Progress in the study on mechanics problems of functionally graded materials and structures [J]. Advances in Mechanics, 2010, 40(5): 528 − 541. (in Chinese)
    [8] Jha D K, Kant T, Singh R K. A critical review of recent research on functionally graded plates [J]. Composite Structures, 2013, 96: 833 − 849. doi:  10.1016/j.compstruct.2012.09.001
    [9] Sankar B V. An elasticity solution for functionally graded beams [J]. Composites Science and Technology, 2001, 61(5): 689 − 696. doi:  10.1016/S0266-3538(01)00007-0
    [10] Zhong Z, Yu T. Analytical solution of a cantilever functionally graded beam [J]. Composites Science and Technology, 2007, 67(3/4): 481 − 488. doi:  10.1016/j.compscitech.2006.08.023
    [11] Li X. A unified approach for analyzing static and dynamic behaviors of functionally graded Timoshenko and Euler–Bernoulli beams [J]. Journal of Sound and Vibration, 2008, 318(4/5): 1210 − 1229. doi:  10.1016/j.jsv.2008.04.056
    [12] Simşek M. Fundamental frequency analysis of functionally graded beams by using different higher-order beam theories [J]. Nuclear Engineering and Design, 2010, 240(4): 697 − 705. doi:  10.1016/j.nucengdes.2009.12.013
    [13] Chakraborty A, Gopalakrishnan S, Reddy J N. A new beam finite element for the analysis of functionally graded materials [J]. International Journal of Mechanical Sciences, 2003, 45(3): 519 − 539. doi:  10.1016/S0020-7403(03)00058-4
    [14] Kadoli R, Akhtar K, Ganesan N. Static analysis of functionally graded beams using higher order shear deformation theory [J]. Applied Mathematical Modelling, 2008, 32(12): 2509 − 2525. doi:  10.1016/j.apm.2007.09.015
    [15] Filippi M, Carrera E, Zenkour A M. Static analyses of FGM beams by various theories and finite elements [J]. Composites Part B: Engineering, 2015, 72: 1 − 9.
    [16] Pradhan S C, Murmu T. Thermo-mechanical vibration of FGM sandwich beam under variable elastic foundations using differential quadrature method [J]. Journal of Sound and Vibration, 2009, 321(1-2): 342 − 362. doi:  10.1016/j.jsv.2008.09.018
    [17] Jiao P F, Wang Y Z, Xu G N, et al. Linear bending of functionally graded beams by differential quadrature method [C]// IOP Conference Series: Earth and Environmental Science. IOP publishing, 2018, 170: 022160.
    [18] 李世荣, 刘平. 功能梯度梁与均匀梁静动态解间的相似转换[J]. 力学与实践, 2010, 32(5): 45 − 49.

    Li Shirong, Liu Ping. Analogous transformation of static and dynamic solutions between functionally graded material beams and uniform beams [J]. Mechanics in Engineering, 2010, 32(5): 45 − 49. (in Chinese)
    [19] 王瑄, 李世荣. 功能梯度Levinson梁自由振动响应的均匀化和经典化表示[J]. 振动与冲击, 2017, 36(18): 70 − 77.

    Wang Xuan, Li Shirong. Homogenized and classical expression for the response of free vibration of simply supported FGM Levinson beams [J]. Journal of Vibration and Shock, 2017, 36(18): 70 − 77. (in Chinese)
    [20] Murín J, Kutiš V. An effective solution of the composite (FGM's) beam structures [J]. Engineering Mechanics, 2008, 15(2): 115 − 132.
    [21] Alshorbagy A E, Eltaher M A, Mahmoud F F. Free vibration characteristics of a functionally graded beam by finite element method [J]. Applied Mathematical Modelling, 2011, 35(1): 412 − 425. doi:  10.1016/j.apm.2010.07.006
    [22] 汪亚运, 彭旭龙, 陈得良. 轴向功能梯度变截面梁的自由振动研究[J]. 固体力学学报, 2015, 36(5): 384 − 391.

    Wang Yayun, Peng Xulong, Chen Deliang. Free vibration of axially functionally graded beams with non-uniform cross-section [J]. Chinese Journal of Solid Mechanics, 2015, 36(5): 384 − 391. (in Chinese)
    [23] Su C, Fan X, Ma H, et al. Green's function method for stability analysis of stochastic structures [J]. Journal of Engineering Mechanics, 2015, 141(3): 4014121. doi:  10.1061/(ASCE)EM.1943-7889.0000842
    [24] Su C, Han D. Elastic analysis of orthotropic plane problems by the spline fictitious boundary element method [J]. Applied Mathematics and Mechanics, 2002, 23(4): 446 − 453.
  • [1] 苏小卒, 王伟.  “依赖应变历史材料”结构动力松弛法静力分析中规避虚假应变历史的非线性弹性增量算法 . 工程力学, doi: 10.6052/j.issn.1000-4750.2019.04.0185
    [2] 何昊南, 于开平.  考虑热对材料参数影响的FGM梁热后屈曲特性研究 . 工程力学, doi: 10.6052/j.issn.1000-4750.2018.03.0142
    [3] 毛丽娟, 马连生.  非均匀热载荷作用下功能梯度梁的非线性静态响应 . 工程力学, doi: 10.6052/j.issn.1000-4750.2015.10.0868
    [4] 丁勇, 黄奇, 黄剑源.  连续桥面简支梁桥静动力特性的理论分析方法研究 . 工程力学, doi: 10.6052/j.issn.1000-4750.2014.01.0096
    [5] 丁 勇, 谢 旭, 黄剑源.  考虑车轮滚动轨迹的桥头跳车动力荷载计算方法及影响因素分析 . 工程力学, doi: 10.6052/j.issn.1000-4750.2011.07.0450
    [6] 朱明亮, 董石麟.  向量式有限元在索穹顶静力分析中的应用 . 工程力学, doi: 10.6052/j.issn.1000-4750.2010.12.0874
    [7] 白兴兰, 黄维平.  深水钢悬链线立管非线性有限元静力分析 . 工程力学,
    [8] 王敏娟, 陈建军.  随机参数智能梁结构闭环控制系统动力响应分析 . 工程力学,
    [9] 项 松, 王克明, 石 宏.  基于逆多元二次径向基函数的复合材料层合板静力分析 . 工程力学,
    [10] 代学灵, 贾宏伟, 李 珠.  某焦煤车间钢梁的耐撞性分析 . 工程力学,
    [11] 王克海, 李 茜.  桥梁抗震的研究进展 . 工程力学,
    [12] 张文福, 刘文洋, 赵文艳, 姚芳.  三向类网架结构弯曲和振动分析的分解刚度法 . 工程力学,
    [13] 王丽, 周锡元, 闫维明.  曲线梁桥地震响应的简化分析方法 . 工程力学,
    [14] 龚耀清, 白俊英.  矩阵位移法分析弹性地基上空间巨型框架结构 . 工程力学,
    [15] 张志宏, 张明山, 董石麟.  平衡矩阵理论的探讨及一索杆梁杂交空间结构的静力和稳定性分析 . 工程力学,
    [16] 肖勇刚, 傅衣铭, 查旭东.  考虑地基耦合效应的矩形中厚板的非线性静力分析 . 工程力学,
    [17] 张楠, 夏禾, DE ROECK Guido.  高速铁路铰接式列车车桥系统动力响应分析 . 工程力学,
    [18] 戴君, 陈建军, 马洪波, 崔明涛.  随机参数结构在随机荷载激励下的动力响应分析 . 工程力学,
    [19] 唐建民, 李长慧, 卓家寿.  拉索穹顶结构几何大变形混合有限元静力分析 . 工程力学,
    [20] 曲华.  贮液构筑物-地基相互作用时的动力特性分析 . 工程力学,
  • 加载中
计量
  • 文章访问数:  40
  • HTML全文浏览量:  28
  • PDF下载量:  2
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-10-31
  • 修回日期:  2020-03-10
  • 网络出版日期:  2020-06-02

功能梯度梁静力与动力分析的格林函数法

doi: 10.6052/j.issn.1000-4750.2019.10.0615
    基金项目:  国家自然科学基金项目(51678252);广州市科学研究计划重点项目(201804020069)
    作者简介:

    邓先琪(1989−),男,广东人,博士生,主要从事结构不确定性问题研究(E-mail: dengxqi@126.com)

    马海涛(1962−),男,辽宁人,教授,博士,博导,主要从事计算力学、结构分析与优化理论与算法研究(E-mail: htma@gzhu.edu.cn)

    通讯作者: 苏 成(1968−),男,广东人,教授,博士,博导,主要从事结构随机振动和计算力学研究(E-mail: cvchsu@scut.edu.cn)
  • 中图分类号: O342

摘要: 功能梯度梁静动态响应的数值分析方法一般局限于有限元法,存在有限元法的固有缺点,有必要发展新的数值求解方法。将功能梯度梁静力分析的控制微分方程转化为与匀质材料梁静力分析控制微分方程相一致的形式,并利用匀质材料梁静力问题的格林函数,开展功能梯度梁的静力分析。在此基础上,进一步推导获得功能梯度梁的柔度矩阵,据此建立功能梯度梁的运动方程,开展功能梯度梁的动力特性分析和动力响应分析。数值算例表明,采用格林函数法可以高效准确地分析功能梯度梁的静力响应与动力行为,验证了方法的计算精度与效率。

English Abstract

邓先琪, 苏成, 马海涛. 功能梯度梁静力与动力分析的格林函数法[J]. 工程力学. doi: 10.6052/j.issn.1000-4750.2019.10.0615
引用本文: 邓先琪, 苏成, 马海涛. 功能梯度梁静力与动力分析的格林函数法[J]. 工程力学. doi: 10.6052/j.issn.1000-4750.2019.10.0615
Xian-qi DENG, Cheng SU, Hai-tao MA. GREEN’S FUNCTION METHOD FOR STATIC AND DYNAMIC ANALYSIS OF FUNCTIONALLY GRADED BEAMS[J]. Engineering Mechanics. doi: 10.6052/j.issn.1000-4750.2019.10.0615
Citation: Xian-qi DENG, Cheng SU, Hai-tao MA. GREEN’S FUNCTION METHOD FOR STATIC AND DYNAMIC ANALYSIS OF FUNCTIONALLY GRADED BEAMS[J]. Engineering Mechanics. doi: 10.6052/j.issn.1000-4750.2019.10.0615
  • 功能梯度材料(Functionally graded material)[1]指由两种或两种以上不同材料复合而成,且各组成成分的结构和性能在空间上沿着某一方向呈连续梯度变化的一种新型非匀质复合材料。由于其新颖的材料可设计性思想和优良的结构性能,近年来引起学者的广泛关注和研究[2-5]。文献[6-8]对功能梯度材料的力学性能与研究进展进行了综述。

    目前,关于材料性能沿梁高度方向梯度变化的功能梯度梁静动态响应分析问题,解析方法和数值方法已有大量的研究。Sankar[9]假定材料弹性模量沿梁高度方向按指数函数变化,采用Euler-Bernoulli梁理论给出了功能梯度简支梁在正弦分布荷载作用下的弹性力学解;而Zhong和Yu[10]则基于Airy应力函数给出了功能梯度悬臂梁的二维弹性力学解。Li[11]提出了功能梯度Timoshenko梁和Euler-Bernoulli梁的静力与动力分析统一列式。Simsek[12]利用不同的高阶梁理论分析了功能梯度梁的基本频率。Chakraborty等[13]基于一阶剪切变形理论提出了一种新的梁单元来研究功能梯度梁的静力、自由振动和波动力学行为;而Kadoli等[14]则基于高阶剪切变形理论,提出了高阶梁单元来进行功能梯度梁的静力分析。Filippi等[15]基于一类统一列式提出了一种功能梯度梁单元模型,并对功能梯度梁静力问题做了相关的研究。Pradhan和Murmu[16]采用微分求积法研究了材料参数沿高度方向呈幂函数形式分布的弹性地基功能梯度梁和层合梁的振动问题。Jiao等[17]基于一阶剪切变形与经典梁理论,采用微分求积法分析了功能梯度梁的线性弯曲问题。李世荣和刘平[18]基于Euler-Bernoulli梁理论,比较了功能梯度梁和均匀梁的控制微分方程,得到了功能梯度梁与均匀梁静动态解之间的相似转换关系。王瑄和李世荣[19]基于Levinson高阶剪切变形理论,研究了两端简支功能梯度梁自由振动频率与经典梁理论下参考匀质梁自由振动频率之间的解析转换关系。以上文献中仅考虑了材料沿梁高度梯度变化的情况。

    相比于材料性能沿梁高度方向梯度变化的功能梯度梁,关于轴向功能梯度梁静力与动力分析的研究相对较少。轴向功能梯度梁的控制微分方程本质上是一个变系数常微分方程,通常情况下很难获得该类方程的解析解。在数值方法方面,Murín和Kutiš[20]采用传递系数描述功能梯度梁在空间上变化的弹性模量,建立了功能梯度Euler-Bernoulli梁的等效刚度矩阵。Alshorbagy等[21]提出了材料沿轴向或高度方向梯度变化的功能梯度梁动力特性分析的有限元法。汪亚运等[22]基于切比雪夫多项式展开技术研究了轴向功能梯度变截面梁的自由振动问题。

    本文主要研究弹性模量沿轴向梯度变化的功能梯度Euler-Bernoulli等截面梁的静力与动力分析问题。由于这种梁的弹性模量沿轴向呈函数变化,因此采用传统有限元法进行结构分析时,不可避免地要进行有限单元网格细分,随之带来的问题是计算规模的增大与计算效率的降低。为了提高功能梯度梁静动力问题的计算精度与效率,本文首先将功能梯度梁静力分析的控制微分方程转化为与匀质材料等截面梁静力分析控制微分方程相一致的形式,然后利用无限长匀质材料等截面梁静力问题的格林函数[23],提出了功能梯度梁静力分析的格林函数法。在此基础上,进一步推导获得功能梯度梁的柔度矩阵,据此建立功能梯度梁的运动方程,开展功能梯度梁的动力特性分析和动力响应分析。

    • 考虑弹性模量沿轴向梯度变化的功能梯度梁静力分析的控制微分方程为:

      $$\left\{ \begin{aligned}& \frac{{\rm{d}}}{{{\rm{d}}x}}\left( {\eta (x)EA\frac{{{\rm{d}}u(x)}}{{{\rm{d}}x}}} \right) = - p(x) \\& \frac{{{{\rm{d}}^2}}}{{{\rm{d}}{x^2}}}\left( {\eta (x)EI\frac{{{{\rm{d}}^2}v(x)}}{{{\rm{d}}{x^2}}}} \right) = q(x) \end{aligned} \right.$$ (1)

      式中:$x$为梁的轴向坐标;$u(x)$$v(x)$分别为梁的轴向位移和横向位移;$E$为参考匀质梁的弹性模量;$\eta (x)$为反映弹性模量沿轴向变化的函数;$A$$I$分别为梁的截面面积和截面惯性矩;$p(x)$$q(x)$分别为梁的轴向和横向分布荷载。

      根据微分运算关系,式(1)可以转化为:

      $$\left\{\begin{align} & EA\frac{{{{\rm{d}}^2}u(x)}}{{{\rm{d}}{x^2}}} = - {F_1}(x) \\ & EI\frac{{{{\rm{d}}^4}v(x)}}{{{\rm{d}}{x^4}}} = {F_2}(x) \\ \end{align} \right.$$ (2)

      式中,${F_1}(x)$${F_2}(x)$分别为梁的等效分布荷载,表示为:

      $$\left\{ \begin{align} {F_1}(x) = & a(x)p(x) + b(x)EA\frac{{{\rm{d}}u(x)}}{{{\rm{d}}x}} \\ {F_2}(x) = & a(x)q(x) - 2b(x)EI\frac{{{{\rm{d}}^3}v(x)}}{{{\rm{d}}{x^3}}} - \\ &c(x)EI\frac{{{{\rm{d}}^2}v(x)}}{{{\rm{d}}{x^2}}} \\ \end{align} \right.$$ (3)

      式中:

      $$\left\{ \begin{aligned}& a(x) = \frac{1}{{\eta (x)}} \\& b(x) = \frac{1}{{\eta (x)}}\frac{{{\rm{d}}\eta {\rm{(}}x{\rm{)}}}}{{{\rm{d}}x}} \\& c(x) = \frac{1}{{\eta (x)}}\frac{{{{\rm{d}}^2}\eta (x)}}{{{\rm{d}}{x^2}}} \end{aligned} \right.$$ (4)

      由式(2)可以看出,功能梯度梁静力分析的控制微分方程可以写成与匀质梁静力分析的控制微分方程相似的形式,唯一的差别在于式(2)中的等效分布荷载与方程左端梁的位移相关联。尽管如此,仍可以利用匀质梁静力问题的格林函数[23]进行求解。

      利用匀质梁静力问题的格林函数进行求解时,在内力计算上存在一些不同。考虑到材料的非匀质性,功能梯度梁内力与位移之间的微分关系为:

      $$\left\{ \begin{align} N(x) = & \eta (x)EA\frac{{{\rm{d}}u(x)}}{{{\rm{d}}x}} \\ M(x) = & \eta (x)EI\frac{{{{\rm{d}}^2}v(x)}}{{{\rm{d}}{x^2}}} \\ Q(x) = & \frac{{\rm{d}}}{{{\rm{d}}x}}\left( {\eta (x)EI\frac{{{{\rm{d}}^2}v(x)}}{{{\rm{d}}{x^2}}}} \right) = \\ & \eta (x)EI\frac{{{{\rm{d}}^3}v(x)}}{{{\rm{d}}{x^3}}} + b(x)M(x) \\ \end{align} \right.$$ (5)

      式中,$b(x)$见式(4)。

      由式(5)可见,在功能梯度梁的静力分析中,内力与位移微分关系在形式上与匀质梁问题的内力与位移微分关系并不完全一致。因此,在采用匀质梁静力问题的格林函数进行求解时,需要考虑函数$\eta (x)$及修正项$b(x)M(x)$对内力的影响。

      利用式(5),可将式(3)给出的等效荷载表示为:

      $$\left\{ \begin{align} {F_1}(x) = & a(x)( {p(x) + b(x)N(x)} ) \\ {F_2}(x) = & a(x)[ {q(x) - 2b(x)Q(x) + } \\ & {( {2{b^2}(x) - c(x)} )M(x)}] \\ \end{align} \right.$$ (6)

      由式(6)可见,等效荷载可以写成与内力$N(x)$$Q(x)$$M(x)$相关的形式,因此在采用格林函数法分析时,需要建立关于内力响应的补充方程,方可求解。

    • 利用功能梯度梁静力分析控制微分方程和匀质梁静力分析控制微分方程在形式上的相似性,可以将功能梯度梁转化为如图1所示的等效匀质梁ab进行求解。其中,$x - y$为梁的坐标系,梁端点轴向坐标为${x_{\rm{a}}}$${x_{\rm{b}}}$,弹性模量取为参考匀质梁的弹性模量$E$,截面面积和截面惯性矩分别取为功能梯度梁的截面面积$A$和截面惯性矩$I$,轴向和横向分布荷载分别取为式(6)所示的等效分布荷载${F_1}(x)$${F_2}(x)$

      图  1  与功能梯度梁等效的匀质梁

      Figure 1.  A uniform beam equivalently to the functionally graded beam

      采用域外虚荷载方法[24],将等效匀质梁ab嵌入到具有相同材料性质与截面性质的无限长梁中,并在梁ab域外${\xi _j}(j = 1,2)$处布设6个大小未知的集中虚荷载$X_j^i\;(i = 1,2,3;j = 1,2)$。根据叠加原理,在域内等效荷载${F_i}(x)\;(i = 1,2)$与域外虚荷载$X_j^i\;(i = 1,2,3;j = 1,$ $2)$的共同作用下,梁ab上任意一点的位移与内力响应为

      $$ \left\{ \begin{array}{l} u(x) = \displaystyle\sum\limits_{i = 1}^3 {\displaystyle\sum\limits_{j = 1}^2 {G_u^i(x;{\xi _j})X_j^i} } + \displaystyle\sum\limits_{i = 1}^2 {\int_{{x_{\rm{a}}}}^{{x_{\rm{b}}}} {G_u^i(x;\xi ){F_i}(\xi ){\rm{d}}\xi } } \\ v(x) = \displaystyle\sum\limits_{i = 1}^3 {\displaystyle\sum\limits_{j = 1}^2 {G_v^i(x;{\xi _j})X_j^i} } + \displaystyle\sum\limits_{i = 1}^2 {\int_{{x_{\rm{a}}}}^{{x_{\rm{b}}}} {G_v^i(x;\xi ){F_i}(\xi ){\rm{d}}\xi } } \\ \theta (x) = \displaystyle\sum\limits_{i = 1}^3 {\displaystyle\sum\limits_{j = 1}^2 {G_\theta ^i(x;{\xi _j})X_j^i} } + \displaystyle\sum\limits_{i = 1}^2 {\int_{{x_{\rm{a}}}}^{{x_{\rm{b}}}} {G_\theta ^i(x;\xi ){F_i}(\xi ){\rm{d}}\xi } } \\ N(x) = \eta (x)\left[ {\displaystyle\sum\limits_{i = 1}^3 {\displaystyle\sum\limits_{j = 1}^2 {G_N^i(x;{\xi _j})X_j^i} } + } \right. \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\displaystyle\sum\limits_{i = 1}^2 {\int_{{x_{\rm{a}}}}^{{x_{\rm{b}}}} {G_N^i(x;\xi ){F_i}(\xi ){\rm{d}}\xi } } } \right] \\ Q(x) = \eta (x)\left[ {\displaystyle\sum\limits_{i = 1}^3 {\displaystyle\sum\limits_{j = 1}^2 {G_Q^i(x;{\xi _j})X_j^i} } + } \right. \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\displaystyle\sum\limits_{i = 1}^2 {\int_{{x_{\rm{a}}}}^{{x_{\rm{b}}}} {G_Q^i(x;\xi ){F_i}(\xi ){\rm{d}}\xi } } } \right] + b(x)M(x) \\ M(x) = \eta (x)\left[ {\displaystyle\sum\limits_{i = 1}^3 {\displaystyle\sum\limits_{j = 1}^2 {G_M^i(x;{\xi _j})X_j^i} } + } \right. \\ \left. {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\displaystyle\sum\limits_{i = 1}^2 {\int_{{x_{\rm{a}}}}^{{x_{\rm{b}}}} {G_M^i(x;\xi ){F_i}(\xi ){\rm{d}}\xi } } } \right] \\[-12pt] \end{array} \right. $$ (7)

      式中,$G_R^i(R = u,v,\theta ,N,Q,M;\;i = 1,2,3)$为匀质梁静力问题的格林函数[23]。值得注意的是,由于采用了匀质梁静力问题的格林函数,根据式(5)所示的内力与位移微分关系,式(7)中的轴力$N(x)$、剪力$Q(x)$和弯矩$M(x)$表达式必须考虑函数$\eta (x)$和剪力修正项$b(x)M(x)$的影响。

      式(7)中的域内积分项与等效荷载${F_1}(x)$${F_2}(x)$有关,而由式(6)可知,${F_1}(x)$${F_2}(x)$除了与真实荷载$p(x)$$q(x)$有关外,还与梁的内力有关。分别采用$l$次和$m$次多项式表示轴力$N(x)$和弯矩$M(x)$的轴向分布,有:

      $$\left\{ \begin{array}{l} N(x) = \displaystyle\sum\limits_{i = 0}^l {{\alpha _i}{x^i}} \\ M(x) = \displaystyle\sum\limits_{i = 0}^m {{\beta _i}{x^i}} \\ \end{array} \right.$$ (8)

      此外,根据剪力与弯矩的微分关系有:

      $$Q(x) = \frac{{{\rm{d}}M(x)}}{{{\rm{d}}x}} = \sum\limits_{i = 1}^m {i{\beta _i}{x^{i - 1}}} $$ (9)

      于是,等效荷载${F_1}(x)$${F_2}(x)$中与内力相关部分可以改用内力系数${\alpha _i}(i = 0,1,2, \cdots ,l)$${\beta _i}(i = 0,1,$ $2, \cdots ,m)$进行表达。

      为了计算式(7)中的域内积分项,可将图1中的梁ab划分成$n$个子段,并将其中每段的中点坐标记为${x_k}\;(k = 1,2, \cdots ,n)$。将式(6)、式(8)和式(9)代入式(7)中的积分项,并采用分段矩形积分公式计算积分,则式(7)可以表达为:

      $$ { {R}}(x) = { {A}}(x){ {X}} + { {B}}(x) { {f}} + { {C}}(x) { {Z}} $$ (10)

      式中,${ {R}}(x)$${ {X}}$${ {f}}$${ {Z}}$分别为响应向量、虚荷载向量、子段荷载向量和内力系数向量:

      $$\left\{ \begin{array}{l} \ { {R}}(x)\ = [u(x)\;\;v(x)\;\;\theta (x)\;\;N(x)\;\;Q(x) \\ \;\;\;\;\;\;\;\;\;\;\;\;\quad M(x){]^{\rm{T}}} \\ \ { {X}}\ = {[X_1^1\;\;X_1^2\;\;\;X_1^3\;\;X_2^1\;\;X_2^2\;\;X_2^3]^{\rm{T}}} \\ \ { {f}}\ = [p({x_1})\;\;q({x_1})\;\;p({x_2})\;\;q({x_2})\;\;...\;\; \\ \;\;\;\;\;\;\;\;\;\;p({x_n})\;\;q({x_n}){]^{\rm{T}}} \\ \ { {Z}}\ = [{\alpha _0}\;{\alpha _1}\;{\alpha _2}\; \cdots \;{\alpha _l}\;{\beta _0}\;{\beta _1}\;{\beta _2}\;\; \cdots\;{\beta _m}{]^{\rm{T}}} \end{array} \right.$$ (11)

      ${ {A}}(x)$${ {B}}(x)$${ {C}}(x)$分别为对应于虚荷载向量、子段荷载向量和内力系数向量的影响矩阵。

      为了求解未知虚荷载向量${ {X}}$及内力系数向量${ {Z}}$,首先需要考虑梁端边界条件。对于平面梁问题,每个梁端可以列出3个边界条件。因此,将梁端坐标${x_{\rm{a}}}$${x_{\rm{b}}}$代入式(10),由梁端边界条件可得与之对应的6个方程,可写成以下形式:

      $${ {R}}_0 = { {A}}_0{ {X}} + { {B}}_0{ {f}} + { {C}}_0{ {Z}} $$ (12)

      式中:$ {{{R}}_0} $含梁端已知位移和内力分量;${{{A}}_0}$${{{B}}_0}$${{{C}}_0}$分别为与之对应的影响矩阵。

      注意到在式(12)中,除了$ {{X}} $中的6个未知虚荷载以外,还有$ {{Z }}$中的($l + m + 2$)个未知内力系数,因此需要补充建立($l + m + 2$)个方程才能求解。为此,首先根据式(8)可将轴力和弯矩表示为以下形式:

      $$\left\{ {\begin{array}{*{20}{c}} {N(x)} \\ {M(x)} \end{array}} \right\} = \left[ {\begin{array}{*{20}{c}} 1&x& \ldots &{{x^l}}&0&0& \ldots &0 \\ 0&0& \ldots &0&1&x& \ldots &{{x^m}} \end{array}} \right]{Z} $$ (13)

      理论上,此式与式(10)在同一点处应给出相同的内力值。基于此,在梁域内选择($l + 1$)个点,利用轴力$N(x)$相等可以得到($l + 1$)个方程;类似的,在梁域内选择($m + 1$)个点,利用弯矩$M(x)$相等可以得到($m + 1$)个方程。最终可得:

      $${H Z} = {\tilde{ A }}{X}+{\tilde{ B }}{f}+ {\tilde{ C }}{Z} $$ (14)

      式中,${{H}}$的各行取决于式(13)系数矩阵在所选择点处的取值,而$\tilde{ {A}}$$\tilde{{ B}}$$\tilde { {C} }$的各行取决于式(10)中对应影响矩阵在所选择点处的取值。

      联立式(12)和式(14),可以解得:

      $$\left\{ {\begin{aligned} & {X} = { {J}_1} { {{R}_0}} + { {J}_2} { f} \\& { Z} = {{L}_1} {{R}_0} + {{L}_2}{ f} \end{aligned} }\right.$$ (15)

      式中:

      $$\left\{ {\begin{aligned}& {{{J}_1} = {{{{A}_0}}^{ - 1}}({{I}} - {{C}_0}{{L}_1})}\\& {{{J}_2} = - {{{{A}_0}}^{ - 1}}({{B}_0} + {{C}_0}{{L}_2})}\\& {{{L}_1} = {{T}^{ - 1}}\tilde{ {A}}{{{{A}_0}}^{ - 1}}}\\& {{{L}_2} = {{T}^{ - 1}}(\tilde{ {B}} - \tilde {{A}}{{{{A}_0}}^{ - 1}}{{B}_0})}\\& {{T} = {H} + \tilde {{A}}{{{{A}_0}}^{ - 1}}{{C}_0} - \tilde {{C}}} \end{aligned}} \right.$$ (16)

      式中,${{I}}$$6 \times 6$的单位矩阵。

      从以上推导可见,由式(15)求得的基本未知向量${{X }}$$ {{Z}} $分别包含6个虚荷载和($l + m + 2$)个内力系数,基本未知量数目与梁的子段划分无关,即与子段数$n$无关。

      $ {{X}} $$ {{Z}} $代入式(10),即可得到功能梯度梁内任意一点的位移和内力响应为:

      $$ {R}(x) = {\delta _1}(x) {{R}_0} + {\delta _2}(x) {f} $$ (17)

      式中:

      $$\left\{ {\begin{array}{*{20}{l}} {[{\delta _1}(x)] = [A(x)][{J_1}] + [C(x)][{L_1}]}\\ {[{\delta _2}(x)] = [B(x)] + [A(x)][{J_2}] + [C(x)][{L_2}]} \end{array}} \right.$$ (18)
    • $ {{{R}}_0} = 0 $,由式(17)可以得到由子段荷载向量$ {{f }}$引起的梁各子段中点${x_k}(k = 1,2, \cdots ,n)$处的位移为:

      $$ {w} = {\delta} {f} $$ (19)

      式中,${ w} = {[u({x_1})\;v({x_1})\;u({x_2})\;v({x_2})\; \cdots \;u({x_n})\;v({x_n})]^{\rm{T}}}$${ f} $见式(11),${{\delta}} $可由${{{\delta}} _2}(x)$前两行元素在${x_k}\;(k = 1,\; 2\; \cdots ,n)$处取值得到。根据柔度矩阵的定义,式(19)中的${{\delta}} $即为功能梯度梁对应于$ {{w}} $${{ f}} $的柔度矩阵。

    • 为了开展功能梯度梁的动力分析,除了外荷载外,还需要考虑功能梯度梁的惯性力和阻尼力。于是,若将质量矩阵和阻尼矩阵记为${{M}}$${D}$,在外荷载$\{ f(t)\} $、惯性力$( - {{M}}\{ \ddot w(t)\} )$和阻尼力$( -{{ D}}\{ \dot w(t)\} )$的共同作用下,由式(19)可得:

      $$\{ w(t)\} = {\delta} (\{ f(t)\} - {M}\{ \ddot w(t)\} - {D}\{ \dot w(t)\} )$$ (20)

      整理后得用柔度矩阵表示的功能梯度梁的运动方程:

      $$ {\delta M}\{ \ddot w(t)\} + {\delta D}\{ \dot w(t)\} + \{ w(t)\} = {\delta }\{ f(t)\} $$ (21)

      将各子段梁的分布质量堆聚在各子段中点${x_k}\;(k = 1,2, \cdots ,n)$上,把各子段堆聚质量记为${m_k}(k = 1,2, \cdots ,n)$,则质量矩阵为:

      $${M }= \left[ {\begin{array}{*{20}{l}} {{m_1}}&{}&{}&{}&{}&{}&{} \\ {}&{{m_1}}&{}&{}&{}&0&{} \\ {}&{}&{{m_2}}&{}&{}&{}&{} \\ {}&{}&{}&{{m_2}}&{}&{}&{} \\ {}&{}&{}&{}& \ddots &{}&{} \\ {}&0&{}&{}&{}&{{m_n}}&{} \\ {}&{}&{}&{}&{}&{}&{{m_n}} \end{array}} \right]$$ (22)

      阻尼矩阵可采用Rayleigh阻尼模型生成。

    • 式(21)可以退化为功能梯度梁无阻尼自由振动方程。利用传统的特征值问题求解方法,可以获得功能梯度梁的自振频率和振型。

      在动力荷载$\{ f(t)\} $的作用下,采用传统的振型分解法或直接积分法求解式(21),可以获得功能梯度梁的位移响应$\{ w(t)\} $、速度响应$\{ \dot w(t)\} $和加速度响应$\{ \ddot w(t)\} $。同时,考虑外荷载、惯性力和阻尼力的作用,由式(17)可以进一步获得功能梯度梁内任意一点处的位移和内力响应为:

      $$\begin{split} {R}(x,t)= &[{\delta _1}(x)] {{R}_0} + [{\delta _2}(x)](\{ f(t)\} - \\ & {M}\{ \ddot w(t)\} - {D}\{ \dot w(t)\} ) \end{split} $$ (23)
    • 本文所采用的格林函数法和有限元法计算程序均采用Python编程语言进行编写,并在桌面电脑(Intel®Core™i7-2620 CPU, 8-GB RAM)上运行。

    • 图2所示的一端固支、另一端简支功能梯度梁受轴向三角形分布荷载$p(x) = 10(1 - {x / l})\;{\rm{kN/m}}$与横向三角形分布荷载$q(x) = 10(1 - {x / l})\;{\rm{kN/m}}$共同作用。梁的截面面积为$A = 1.0 \times {10^{ - 3}}\;{{\rm{m}}^2}$,截面惯性矩为$I = 2.083 \times {10^{ - 7}}\;{{\rm{m}}^4}$,长度为$l = 1\;{\rm{m}}$。假定梁的弹性模量沿轴向按以下规律变化:

      $$\begin{split} {E_{\rm{b}}}(x) =& E\eta (x) = 37.7\left( {1 + 0.8284\frac{x}{l} - } \right. \\& \left. {0.3682\frac{{{x^2}}}{{{l^2}}}} \right)\;{\rm{GPa,}}\;\;\;0 \leqslant x \leqslant l \end{split} $$ (24)

      式中,$E$为参考匀质梁的弹性模量,取为$E = 37.7\;{\rm{GPa}}$

      图  2  一端固支、一端简支功能梯度梁

      Figure 2.  A clamped-simply supported functionally graded beam

      分别采用解析法、格林函数法和有限元法计算功能梯度梁的位移和内力。根据材料力学理论,该功能梯度梁的轴向与横向位移、轴力与弯矩解析解表达式分别为:

      $$\left\{ \begin{split} u(x) = &[ - 4.500\ln ( - {x^2} + 2.250x + 2.720) + \\ &65.531{\tanh ^{ - 1}}(0.500x - 0.560) - \\ &36.020x + 46.330] \times {10^{ - 5}}\;\;{\rm{m}} \\ v(x) = &[980.906(x - 1.125) \cdot \\ &\ln ( - {x^2} + 2.250x + 2.715) - \\ &638.230\ln (3.120 - x)- \\ &178.036\ln (x + 0.870) + \\ &409.085(x - {\rm{9}}{\rm{.56}}9) \cdot \\ &{\tanh ^{ - 1}}(0.501x - 0.564) + \\ &960.690{x^3} - 216.185{x^2} - 3203.006x + \\ &4302.863] \times {10^{ - 3}}\;\;{\rm{m}} \\ \end{split} \right.$$ (25)

      $$\left\{ \begin{split} N(x) = & 5{x^2} - 10x + 3.166\;\;{\rm{kN}} \\ M(x) = & - 1.667{x^3} + 5{x^2} - 3.959x +\\ & 0.625\;\;{\rm{kN}} \cdot {\rm{m}} \end{split} \right.$$ (26)

      在格林函数法中,采用分段矩形积分公式计算式(7)中的域内积分项,分别取$n = 4$$10$$20$$50$$100$以考虑子段数对计算精度的影响。在轴向和横向三角形分布荷载作用下,梁的轴力和弯矩分别为关于$x$的2次和3次多项式,因此分别将式(8)中的轴力和弯矩表示为相应次数的多项式,即$l = 2$$m = 3$,并在梁内均匀选取($l + 1$)和($m + 1$)个点分别建立关于轴力和弯矩的补充方程。将两端虚荷载施加位置与梁左右端点的距离表示为d1d2,分别考虑${d_1} = {d_2} = 0.2\;{\rm{m}}$${d_1} = {d_2} = 0.4\;{\rm{m}}$${d_1} = 0.2\;{\rm{m}}$${d_2} = 0.4\;{\rm{m}}$三种情况对格林函数法计算结果的影响。

      在有限元法中,将功能梯度梁离散成${N_{\rm e}}$个基于Hermite插值的平面二节点梁单元,分别取${N_{\rm e}} = 4$$10$$20$$50$$100$以考虑离散单元数对计算精度的影响,其中每个梁单元的弹性模量取为梁单元中点处的弹性模量值。

      表1给出了功能梯度梁跨中点C的轴向位移和横向位移结果,其收敛情况如图3所示。表2给出了功能梯度梁固定端点A的轴力和弯矩结果,其收敛情况如图4所示。由上述图表可见,当子段数$n \geqslant 20$时,格林函数法的计算结果与解析解高度吻合,说明所提出方法具有理想的计算精度。此外,由计算结果还可以看出,虚荷载位置的三种选择对格林函数法计算结果没有任何影响,说明所提出方法具有良好的数值稳定性。

      表 1  功能梯度梁跨中点C的轴向与横向位移

      Table 1.  Axial and transverse displacement of node C at midspan of the functionally graded beam

      位移方法子段数n,离散单元数Ne
      4102050100
      轴向位移/
      (×10−5 m)
      解析解 1.328
      有限元法 1.314 1.326 1.327 1.328 1.328
      格林函数法
      (d1=d2=0.2 m)
      1.349 1.331 1.329 1.328 1.328
      格林函数法
      (d1=d2=0.4 m)
      1.349 1.331 1.329 1.328 1.328
      格林函数法
      (d1=0.2 m, d2=0.2 m)
      1.349 1.331 1.329 1.328 1.328
      横向位移/
      (×10−3 m)
      解析解 −2.447
      有限元法 −2.403 −2.439 −2.445 −2.447 −2.447
      格林函数法
      (d1=d2=0.2 m)
      −2.489 −2.451 −2.448 −2.447 −2.447
      格林函数法
      (d1=d2=0.4 m)
      −2.489 −2.451 −2.448 −2.447 −2.447
      格林函数法
      (d1=0.2 m, d2=0.2 m)
      −2.489 −2.451 −2.448 −2.447 −2.447

      表 2  功能梯度梁固定端点A的轴力与弯矩

      Table 2.  Axial force and bending moment of node A at fixed end of the functionally graded beam

      内力方法子段数n,离散单元数Ne
      4102050100
      轴力/kN 解析解 3.166
      有限元法 3.181 3.168 3.166 3.166 3.166
      格林函数法
      (d1=d2=0.2 m)
      3.066 3.149 3.161 3.165 3.165
      格林函数法
      (d1=d2=0.4 m)
      3.066 3.149 3.161 3.165 3.165
      格林函数法
      (d1=0.2 m, d2=0.2 m)
      3.066 3.149 3.161 3.165 3.165
      弯矩/(kN$ \cdot $m) 解析解 −0.625
      有限元法 −0.634 −0.627 −0.626 −0.625 −0.625
      格林函数法
      (d1=d2=0.2 m)
      −0.677 −0.634 −0.628 −0.626 −0.625
      格林函数法
      (d1=d2=0.4 m)
      −0.677 −0.634 −0.628 −0.626 −0.625
      格林函数法
      (d1=0.2 m, d2=0.2 m)
      −0.677 −0.634 −0.628 −0.626 −0.625

      图  3  功能梯度梁跨中点C的轴向与横向位移

      Figure 3.  Axial and transverse displacement of node C at midspan of the functionally graded beam

      图  4  功能梯度梁固定端点A的轴力和弯矩

      Figure 4.  Axial force and bending moment of node A at fixed end of the functionally graded beam

      在计算效率方面,格林函数法中的基本未知量为虚荷载$\{ X\} $和内力系数$\{ Z\} $,详见式(12)和式(14),未知量数目并不受子段数$n$的影响,因此即使$n$取较大值也不会影响方程组的求解效率。而在有限元法中,随着离散单元数${N_{\rm e}}$的增加,基本未知量数目和刚度矩阵规模也相应增多,从而计算工作量也会显著增加,计算时间延长。

    • 考虑图2所示一端固支、另一端简支的功能梯度梁。梁的截面面积、截面惯性矩和长度与算例3.1中相应参数相同。假设梁的弹性模量和密度均沿轴向呈3次多项式函数形式变化,表达式为:

      $$\left\{ \begin{aligned}& {E_{\rm{b}}}(x) = ({E_{\rm{L}}} - {E_{\rm{R}}}){\left( {1 - \frac{x}{l}} \right)^3} + {E_{\rm{R}}} \\& {\rho _{\rm{b}}}(x) = ({\rho _{\rm{L}}} - {\rho _{\rm{R}}}){\left( {1 - \frac{x}{l}} \right)^3} + {\rho _{\rm{R}}} \end{aligned} \right.,\;\;0 \leqslant x \leqslant l$$ (27)

      式中:${E_{\rm{L}}} = 210\;{\rm{GPa}}$${E_{\rm{R}}} = 390\;{\rm{GPa}}$分别表示梁左、右两端处弹性模量的取值;${\rho _{\rm{L}}} = 7800\;{\rm{kg/}}{{\rm{m}}^3}$${\rho _{\rm{R}}} = 3960\;{\rm{kg/}}{{\rm{m}}^3}$分别表示梁左、右两端处密度的取值。

      分别采用格林函数法和有限元法计算功能梯度梁的自振频率。在格林函数法中,反映虚荷载位置的距离取为${d_1} = {d_2} = 0.2\;{\rm{m}}$。为了考察子段数$n$对计算精度的影响,分别取$n = 10$$20$$30$$40$$50$$200$。由于梁在自由振动过程中其内力分布不再满足简单的多项式分布形式,为了考察弯矩多项式次数$m$对横向自由振动固有频率的影响,分别取$m = 4$、6、8和10,而轴力多项式次数则取为$l = 10$。在有限元法中,将功能梯度梁离散成${N_{\rm{e}}}$个基于Hermite插值的平面二节点梁单元,采用一致质量矩阵。分别取${N_{\rm{e}}} = 10$$20$$30$$40$$50$$200$以考察离散单元数${N_{\rm{e}}}$对计算精度的影响,其中每个梁单元的弹性模量取为梁单元中点处的弹性模量值。

      表3图5给出了功能梯度梁的前三阶横向自由振动固有频率的计算结果。由计算结果可见,当$m > 6$时,在子段数$n$和离散单元数${N_{\rm{e}}}$取为相同的情况下,格林函数法计算得到的前三阶横向自由振动固有频率总是比有限元法计算得到的结果精度更高,收敛速度更快。也就是说,当弯矩多项式次数足以反映真实弯矩分布时,要达到相同的精度,格林函数法所需要的子段数$n$少于有限元法所需要的离散单元数${N_{\rm{e}}}$,因此其计算效率更高。

      表 3  功能梯度梁前三阶横向自由振动固有频率

      Table 3.  The first three natural frequencies of functionally graded beam’s transverse free vibration

      频率方法子段数n,离散单元数Ne
      1020304050200
      第一阶频率/
      (×103 rad/s)
      有限元法1.3201.3231.3241.3241.3241.324
      格林函数法
      (m=4)
      1.3221.3231.3231.3231.3231.323
      格林函数法
      (m=6)
      1.3231.3241.3241.3241.3241.324
      格林函数法
      (m=8)
      1.3231.3241.3241.3241.3241.324
      格林函数法
      (m=10)
      1.3221.3241.3241.3241.3241.324
      第二阶频率/
      (×103 rad/s)
      有限元法4.2754.3074.3124.3144.3154.317
      格林函数法
      (m=4)
      4.2984.3034.3034.3044.3044.304
      格林函数法
      (m=6)
      4.3084.3144.3154.3154.3164.316
      格林函数法
      (m=8)
      4.3084.3154.3164.3174.3174.317
      格林函数法
      (m=10)
      4.3094.3154.3164.3174.3174.317
      第三阶频率/
      (×103 rad/s)
      有限元法8.8538.9929.0189.0279.0319.038
      格林函数法
      (m=4)
      8.8858.8978.8988.8998.8998.899
      格林函数法
      (m=6)
      8.9898.9958.9948.9948.9948.993
      格林函数法
      (m=8)
      9.0119.0329.0349.0359.0359.036
      格林函数法
      (m=10)
      9.0109.0339.0369.0379.0379.038

      图  5  功能梯度梁前三阶横向自由振动固有频率

      Figure 5.  The first three natural frequencies of functionally graded beam’s transverse free vibration

    • 考虑图2所示一端固支、另一端简支的功能梯度梁,其几何参数和材料参数与算例3.2中的相应参数相同。梁受横向均布突加荷载$q(t)$作用,荷载随时间变化曲线如图6所示。选取Rayleigh阻尼模型,阻尼比取为$\zeta = 0.05$(对应的自振频率${\omega _1} = 1324.42\;{\rm{rad/s}}$${\omega _2} = 4317.06\;{\rm{rad/s}}$)。

      图  6  荷载-时间曲线

      Figure 6.  Load-time curve

      分别采用格林函数法和有限元法计算功能梯度梁的位移、速度和加速度。在格林函数法中,反映虚荷载位置的距离取为${d_1} = {d_2} = 0.2\;{\rm{m}}$,内力多项式次数取为$l = m = 8$。为了考察子段数$n$对计算精度的影响,分别取$n = 10$$20$$30$。在有限元法中,将功能梯度梁离散成${N_{\rm{e}}}$个基于Hermite插值的平面二节点梁单元,采用一致质量矩阵。取${N_{\rm{e}}} = 200$以获得相对精确解,其中每个梁单元的弹性模量取为梁单元中点处的弹性模量值。两种方法均采用Newmark-$\beta $积分格式求解运动方程,积分步长均取为$\Delta t = 1.0 \times {10^{ - 4}}\;{\rm{s}}$,共计算$500$个时间步。

      图7给出了功能梯度梁跨中点C的横向位移、速度和加速度结果,上述结果在0.03 s后的时间内基本不再变化。从图中可见,当子段数$n = 30$时,格林函数法计算得到的跨中点C的横向位移、速度和加速度响应与有限元法结果高度吻合,再次说明所提出方法具有理想的计算精度。

      图  7  功能梯度梁跨中点C横向位移、速度和加速度

      Figure 7.  Transverse displacement, velocity, acceleration of node C at midspan of the functionally graded beam

    • 本文提出了一种轴向功能梯度材料Euler-Bernoulli梁静力与动力分析的格林函数法。将功能梯度梁静力分析的控制微分方程转化为与匀质材料梁静力分析控制微分方程相一致的形式,从而可以利用匀质材料梁静力问题的格林函数进行求解。在此基础上,进一步可以推导得到功能梯度梁的柔度矩阵,并建立功能梯度梁的运动方程,开展功能梯度梁动力特性分析与动力响应分析。数值算例表明,对于功能梯度梁的静力和动力问题,格林函数法比有限元法具有更高的计算精度和计算效率。

目录

    /

    返回文章
    返回