引 言
纳米科技被称为21世纪三大尖端技术之一，近几十年来得到迅速发展．纳米材料作为纳米科技研究的主要对象，指的是颗粒或尺寸至少在一维尺度上小于10 nm，具有不同于块状材料性能的一类材料，其中纳米材料的力学性能在其应用中占据了重要位置，吸引国内外学者开展了大量研究，这也促使各种修正的连续介质力学理论得以发展和完善，包括微极与微态理论[1]、偶应力理论[2]、应变梯度理论[3]、非局部理论[4]等．比如，徐晓建和邓子辰[5]利用非局部理论并考虑微结构表面效应，研究了表面吸附物对微纳传感器振动特性的影响．2015年，Lim等[6]提出了非局部应变梯度理论，不同于非局部理论仅考虑应变自身的非局部效应，以及应变梯度理论仅考虑应变梯度的局部效应，该理论基于新的核函数同时计及了应变和应变梯度的非局部效应．纳米材料中波传播的研究表明，基于非局部应变梯度理论的结果与晶格动力学和实验结果更吻合[6]．因此，自该理论提出以来，被大量应用于纳米材料和结构的研究中[7-15]．
随着科技的进步，单一纳米材料的性能已不能完全满足某些复杂或极端环境下的应用要求．为解决这一类问题，研究者们引入了功能梯度的概念，功能梯度纳米材料被成功设计和制造．因其优异的物理、化学、力学、热学等性能，功能梯度纳米材料在多个领域应用广泛，比如纳米机器人[16]、航天器防护结构[17]等．通过各种高阶连续介质力学理论研究功能梯度纳米材料的力学性能，包括弯曲、屈曲、波传播、振动等方面[7-8,12-13,15,18-20]，是近年来的研究热点．比如，Li等[18]基于非局部应变梯度理论和Euler-Bernoulli梁理论，通过广义差分法推导了轴向功能梯度纳米梁的弯曲、屈曲和振动方程并研究尺度依赖特性．Mahinzare等[19]基于非局部应变梯度理论，研究了旋转功能梯度压电纳米圆板在热载荷作用下的振动．Karami和Janghorban[20]研究了由各向异性材料制成的功能梯度纳米梁的自由振动．
功能梯度纳米板结构在工程中的应用潜力巨大，以纳米机器人为例，它是一种根据分子生物学原理设计制造、可在纳米空间进行操作的功能分子器件，在医疗等领域作用重大．纳米机器人中关键部件采用性能更均衡的功能梯度纳米结构，以保证不同环境下的工作可靠性已成为发展趋势[16,21-22]，所以对这类结构的研究具有重要意义．当前，功能梯度纳米材料中关于纳米梁、纳米杆、纳米管等结构的研究较多[7-8,12,15,18,20]，尽管也有一些文献研究功能梯度纳米板或壳结构[13,19]，多采用非局部理论[23-25]，但理论研究功能梯度纳米板弯曲和屈曲，特别是同时揭示挠度与临界载荷等的非局部效应和应变梯度效应的研究还比较缺乏．本文基于近年来提出的非局部应变梯度理论，将纳米板结构与功能梯度材料特性相结合，研究功能梯度纳米板的弯曲和屈曲，为功能梯度纳米板在纳米机器人、航天器等工程中的应用提供了理论基础．
1 非局部应变梯度理论本构
结合传统非局部理论和应变梯度理论而成的非局部应变梯度理论，其积分本构表达如下[6]：
(1)
其中，σ和σ(1)分别表示x处的非局部和高阶非局部应力场，t表示x处的非局部应变梯度全应力场，ε′表示x′处的应变场，C是弹性张量，α0和α1分别是与应变和一阶应变梯度有关的非局部核函数，e0和e1分别是非局部和高阶非局部材料常数，a是非局部特征尺度，l是与高阶应变梯度有关的材料特征尺度，
为梯度算子．因为式(1)中存在关于应变场和应变梯度场的空间积分，该积分本构求解比较困难．与非局部理论积分本构类似，式(1)可在一定条件下转换为如下微分本构[6]：
[1-(ea)22]tij=Cijkl(1-l22)εkl，
(2)
其中，2为Laplace算子，且假设e0=e1=e．
当e=l=0时，式(2)退化为经典理论(CT)应力-应变关系:
tij=Cijklεkl;
(3)
当l=0时，式(2)退化为传统非局部理论(TNT)应力-应变关系：
[1-(ea)22]tij=Cijklεkl;
(4)
当e=0时，式(2)则退化为应变梯度理论(SGT)应力-应变关系：
tij=Cijkl(1-l22)εkl．
(5)
2 物理模型与平衡方程
2.1 功能梯度纳米板的弯曲
图1所示为功能梯度矩形纳米板，长、宽、高分别为La，Lb，h，建立图示坐标系．假设板的上下表面分别为陶瓷和金属，材料的弹性模量E(z)和密度ρ(z)沿厚度方向按幂指数连续变化，则有
(6)
其中，-h/2≤z≤h/2，k为梯度指数，Ec，Em和ρc，ρm分别是陶瓷和金属的弹性模量和密度，其值为：Ec=390 GPa，ρc=3 960 kg/m3，Em=210 GPa，ρm=7 800 kg/m3．
图1 功能梯度纳米板示意图
Fig. 1 Sketch of a functionally graded nanoplate
功能梯度纳米板有一个面内的应力和应变为零，称之为物理中面(z=z0)，具体的定义如下[26]:
(7)
根据经典薄板小挠度理论中的几何方程，并考虑物理中面的概念，得到功能梯度纳米板的应变表达式：
(8)
其中，w为弯曲挠度．结合式(2)，并考虑物理中面的弯矩可由应力沿板厚积分而来，可得
(9)
其中，D=
B(z-z0)2dz是功能梯度纳米板的等效弯曲刚度，
等为弯矩分量，μ为Poisson比．弯曲和屈曲方程可由静力分析推导，也可从更一般的形式中退化得到．比如通过Hamilton原理：
(δUe+δUk-δUs)dt=0，
(10)
其中，δUe，δUk和δUs分别是外力做功、动能和应变能的变分．
应变能变分为
(11)
动能变分为
(12)
其中 I1=ρ(z)dz, I2=(z-z0)2ρ(z)dz．
外力做功变分为
(13)
其中,q为分布载荷．将式(11)～(13)代入式(10)，由δw≠0可得自由振动控制方程：
(14)
将式(9)代入式(14)，得到用弯曲挠度w表达的自由振动控制方程：
(15)
作为特例，弯曲问题中关于时间的导数都为零，可得平衡方程为
2]q．
(16)
分布载荷q可展开为以下形式：
(17)
其中
对于均布载荷，有
(18)
其中，q0代表均布载荷集度，m,n=1,3,5,…．
以四边简支约束边界为例，根据Navier法，将挠度w表达为满足边界条件的级数形式:
(19)
将式(18)和(19)代入式(16)中，可得Cmn,进而弯曲挠度为
(20)
特别地，当ea=l时，非局部应变梯度理论挠度退化为经典理论结果．以往通过数值方法求解纳米梁振动特性时也有类似结论[11]，这是一个颇有趣但有疑惑的结论，因为根据非局部应变梯度本构思想，只有当ea=l=0，即两大内特征参数均为零时，非局部应变梯度理论回归为经典理论．然而，此处解析结果直接显示当ea=l即可退化为经典结果，原因将在后文结合数值算例分析．
对于受均布载荷作用的矩形纳米板，最大挠度在x=La/2、y=Lb/2处，为
(21)
2.2 功能梯度纳米板的屈曲
考虑四边简支功能梯度纳米板，一组对边受有均布压力，在板边的单位长度上的纵向压力为Fx，如图2所示．
图2 四边简支功能梯度纳米板的屈曲受力图
Fig. 2 Buckling diagram of the functionally graded nanoplate with 4 simply supported sides
仍可利用Hamilton原理，此时功能梯度纳米板的外力做功变分表达为
(22)
将式(11)、(12)、(22)代入式(10)中，由δw≠0得屈曲平衡方程，再将式(9)代入其中可得
(23)
对于四边简支边界，将挠度(19)代入式(23)得到
(24)
对于屈曲问题，在满足式(24)的所有纵向载荷中，数值最小的Fx即为临界屈曲载荷．对于给定的材料和结构尺寸，D，La，Lb，l，ea都是确定的．由式(24)可见，若给定m，当n增大时，Fx也随之增大，因此令n=1，可得临界屈曲载荷
(25)
其中
下面讨论m的取值对临界屈曲载荷的影响，因为l,ea在纳米量级，m的取值对λ的影响忽略不计，主要分析m的取值对f的影响．依次令m=1,2,3,…，计算出La/Lb变化时f的数值，如图3所示．
图3 m与La/Lb的值对f的影响
Fig. 3 Influences of m and La/Lb on values of f
由图3可见，邻近两条曲线的交点横坐标分别是
即,当
时，临界屈曲载荷对应m=1；当
时，对应m=2；当
时，对应m=3;……．例如，对功能梯度正方形纳米板La/Lb=1，临界屈曲载荷对应m=1，为
(26)
类似地，当ea=l时，临界屈曲载荷退化为经典理论结果．
3 数值算例与讨论
3.1 弯曲算例与讨论
对于四边简支功能梯度纳米方板，均布载荷取q0=5×107 N/m2．在式(21)中令m，n取到5，事实上更大的取值对最大挠度的影响很小．算例中陶瓷和金属的弹性模量及密度取值同上，Poisson比取μ=0.3，不同理论下的内特征参数取值如表1所示．
令方形纳米板的尺寸为La=Lb=10 nm,h=1 nm，得到4种理论下的最大挠度如图4所示．当梯度指数增加时，各理论下的最大挠度都在增加，且在开始阶段增速较快，之后渐趋缓于一个极值．传统非局部理论下的挠度最高，应变梯度理论下的挠度最低，二者分列经典理论结果的上下侧，可见非局部效应削弱纳米结构等效刚度，而应变梯度效应则增强等效刚度．
表1 4种理论中的内特征参数取值
Table 1 Values of the internal characteristic parameters in 4 theories
theoryea/nml/nmclassical theory00nonlocal theory10strain gradient theory02nonlocal strain gradient theory12
为了分析纳米板尺寸对最大挠度的影响，取ea=l=1 nm；然后分别取La=Lb=10 nm，h=1 nm；La=Lb=10 nm，h=2 nm；La=15 nm，Lb=10 nm，h=1 nm；La=15 nm，Lb=10 nm，h=2 nm；La=Lb=15 nm，h=2 nm，5种情况下结果如图5所示．可以看出：当板厚h与中面的较小尺寸Lb的比值越小(即板越薄)时，最大挠度越大；而板的长度La与宽度Lb相等时，最大挠度相对较小．因此，要提高结构刚度，可考虑使用正方形板且在不影响其他需求情况下增加板厚．
图4 不同理论下的最大挠度对比情况 图5 不同尺寸下纳米板的最大挠度变化情况
Fig. 4 Comparisons of the maximum deflections
Fig. 5 Variations of the maximum deflection under different theories with different sizes
尽管对于给定的材料和结构，非局部尺度参数和材料特征尺度参数都是确定的．但为了研究二者对挠度的影响，令La=Lb=10 nm，h=1 nm，k=1，分别取l=0.5 nm，1 nm，1.5 nm，2 nm，得到最大挠度随非局部尺度参数的变化如图6(a)～(d)所示．可见，在给定的材料特征尺度参数下，最大挠度随着非局部尺度参数的增大而增大，且增速逐渐提高．因此， 非局部效应增大纳米板的弯曲挠度， 体现纳米结构软化行为．分别取ea=0.5 nm，1 nm，1.5 nm，2 nm，计算结果如图7(a)～(d)所示．可见，最大挠度随着材料特征尺度的增大而减小，且减速逐渐变大．因此应变梯度效应降低弯曲挠度，体现硬化行为．综上，非局部应变梯度模型中同时存在软化和硬化现象．
(a) l=0.5 nm (b) l=1 nm
(c) l=1.5 nm (d) l=2 nm
图6 非局部尺度对最大挠度的影响
Fig. 6 Effects of the nonlocal scale on the maximum deflection
3.2 屈曲算例与讨论
首先，令La=Lb=10 nm，h=1 nm，ea=1 nm，分别取l=0 nm，0.5 nm，1 nm，1.5 nm，2 nm，得到临界屈曲载荷与梯度指数的关系如图8所示．随着梯度指数的增大，临界屈曲载荷减小，且开始阶段下降幅度较大，然后逐渐变缓，因此梯度指数增大使得纳米板的力学性能有所下降．随着材料特征尺度的增大，临界屈曲载荷提高．
(a) ea=0.5 nm (b) ea=1 nm
(c) ea=1.5 nm (d) ea=2 nm
图7 材料特征尺度对最大挠度的影响
Fig. 7 Effects of the material characteristic scale on the maximum deflection
图8 不同材料特征尺度参数下临界屈曲载荷随梯度指数的变化图9 临界屈曲载荷随纳米板尺寸的变化
Fig. 8 The critical buckling load vs. the gradient index with Fig. 9 Variations of the critical buckling load different material characteristic scale parameters with respect to sizes of nanoplates
图10 非局部尺度参数对临界屈曲载荷的影响 图11 材料特征尺度参数对临界屈曲载荷的影响
Fig. 10 Influence of the nonlocal scale parameter Fig. 11 Influence of the material characteristic scale on the critical buckling load parameter on the critical buckling load
接着，探讨不同尺寸对临界屈曲载荷的影响，令k=1，ea=1 nm，l=1 nm，分别取La=Lb=10 nm；La=15 nm，Lb=10 nm；La=Lb=15 nm；La=20 nm，Lb=15 nm；La=Lb=20 nm．图9显示临界屈曲载荷随板厚的增加而增大，当然对于薄板来说这个效果是有限的；当长宽比增大时，临界屈曲载荷增大，说明增大长宽比可增强功能梯度纳米板的稳定性．
最后分析非局部和材料特征尺度参数对临界屈曲载荷的影响，令La=Lb=10 nm，h=1 nm，k=1，分别取l=0 nm，0.5 nm，1 nm，1.5 nm，2 nm，计算临界屈曲载荷如图10所示．可以看出，临界屈曲载荷随非局部尺度参数的增大而降低，且降速愈来愈快．分别取ea=0 nm，0.5 nm，1 nm，1.5 nm，2 nm，计算结果如图11所示．可见材料特征尺度参数越大，临界屈曲载荷越大，且增幅越来越大，应变梯度效应增大了临界屈曲载荷．在功能梯度纳米板的临界屈曲载荷中，非局部与应变梯度效应分别导致了结构软化与硬化行为．
图12 非局部与材料特征尺度参数之比与临界屈曲载荷的关系
Fig. 12 Relationship between the ratio of nonlocal to material characteristic scale parameters and the critical buckling load
为了进一步综合衡量非局部与材料特征尺度参数对临界屈曲载荷的影响， 分别令l=0.2 nm， 0.5 nm， 0.8 nm，1 nm，计算临界屈曲载荷随内特征参数之比ea/l的变化，如图12所示．可见，当ea<l时，相同比值ea/l下，增大材料特征尺度参数将提高临界屈曲载荷，此时应变梯度效应占主导作用；当ea>l时，增大材料特征尺度参数将降低临界屈曲载荷，此时非局部效应占主导作用；当ea=l时，正如之前从解析结果推得的结论，无论两大内特征参数如何取值，临界屈曲载荷不变，与经典解一致．非局部尺度与材料特征尺度对纳米结构分别起到了软化和硬化作用，当ea=l时经典结果恢复，证明了非局部软化与应变梯度硬化作用效果在量级上恰好相当，当ea=l时软化与硬化作用相互抵消．当非局部尺度参数与材料特征尺度参数不相等时，尺度数值较大的那个在软/硬化机制中起主导作用．
4 结 论
基于非局部应变梯度理论，结合Kirchhoff薄板模型与物理中面概念，研究了功能梯度纳米板的弯曲与屈曲．通过Navier法得到了挠度和临界屈曲载荷的显式表达，讨论了非局部尺度与材料特征尺度参数、梯度指数、纳米板尺寸对最大挠度和临界屈曲载荷的影响．结论如下：
1) 不同高阶连续介质力学理论下，非局部理论下的挠度相对最大，应变梯度理论下相对最小，非局部应变梯度理论结果介于二者之间．基于非局部应变梯度理论，当功能梯度纳米板厚度增加时，最大挠度减小，临界屈曲载荷增大．
2) 随着梯度指数增大，挠度提高而临界屈曲载荷降低．特别地，当梯度指数在零附近的小范围内变化时，梯度指数对力学性能的影响显著．正方形纳米板的挠度相对较小，而当纳米板长宽比增大时，临界屈曲载荷增大．
3) 最大挠度/临界屈曲载荷随非局部尺度参数的增大而增大/减小， 随材料特征尺度参数的增大而减小/增大．非局部和材料特征尺度参数存在耦合作用， 二者相比数值更大的那个参数在力学机制中占主导作用．
致谢 本文作者衷心感谢暨南大学重大工程灾害与控制教育部重点实验室开放基金(20180930002)对本文的资助．
[1] KAFADAR C B, ERINGEN A C. Micropolar media Ⅰ: the classical theory[J]. International Journal of Engineering Science, 1971, 9(3): 271-305.
[2] TOUPIN R. Elastic materials with couple-stresses[J]. Archive for Rational Mechanics and Analysis, 1962, 11(1): 385-414.
[3] MINDLIN R D, ESHEL N N. On first strain-gradient theories in linear elasticity[J]. International Journal of Solids and Structures, 1968, 4(1): 109-124.
[4] ERINGEN A C, EDELEN D G B. On nonlocal elasticity[J]. International Journal of Engineering Science, 1972, 10(3): 233-248.
[5] 徐晓建, 邓子辰. 非局部因子和表面效应对微纳米材料振动特性的影响[J]. 应用数学和力学, 2013, 34(1): 10-17. (XU Xiaojian, DENG Zichen. Surface effects of adsorption-induced resonance analysis of micro/nanobeams via nonlocal elasticity[J]. Applied Mathematics and Mechanics, 2013, 34(1): 10-17. (in Chinese))
[6] LIM C W, ZHANG G, REDDY J N. A higher-order nonlocal elasticity and strain gradient theory and its applications in wave propagation[J]. Journal of the Mechanics and Physics of Solids, 2015, 78: 298-313.
[7] LI L, LI X, HU Y. Free vibration analysis of nonlocal strain gradient beams made of functionally graded material[J]. International Journal of Engineering Science, 2016, 102: 77-92.
[8] EBRAHIMI F, BARATI M R. A nonlocal strain gradient refined beam model for buckling analysis of size-dependent shear-deformable curved FG nanobeams[J]. Composite Structures, 2017, 159: 174-182.
[9] XU X J, ZHENG M L, WANG X C. On vibrations of nonlocal rods: boundary conditions, exact solutions and their asymptotics[J]. International Journal of Engineering Science, 2017, 119: 217-231.
[10] SAHMANI S, AGHDAM M M. Nonlocal strain gradient shell model for axial buckling and postbuckling analysis of magneto-electro-elastic composite nanoshells[J]. Composites Part B: Engineering, 2018, 132: 258-274.
[11] WANG J, SHEN H M, ZHANG B, et al. Complex modal analysis of transverse free vibrations for axially moving nanobeams based on the nonlocal strain gradient theory[J]. Physica E: Low-Dimensional Systems and Nanostructures, 2018, 101: 85-93.
[12] SAHMANI S, AGHDAM M M, RABCZUK T. Nonlinear bending of functionally graded porous micro/nano-beams reinforced with graphene platelets based upon nonlocal strain gradient theory[J]. Composite Structures, 2018, 186: 68-78.
[13] LU L, ZHU L, GUO X M, et al. A nonlocal strain gradient shell model incorporating surface effects for vibration analysis of functionally graded cylindrical nanoshells[J]. Applied Mathematics and Mechanics(English Edition), 2019, 40(12): 1695-1722.
[14] MIR M, TAHANI M. Graphene-based mass sensors: chaotic dynamics analysis using the nonlocal strain gradient model[J]. Applied Mathematical Modelling, 2020, 81: 799-817.
[15] BARRETTA R, FAGHIDIAN S A, DE SCIARRA F M, et al. On torsion of nonlocal Lam strain gradient FG elastic beams[J]. Composite Structures, 2020, 233: 111550.
[16] LI J X, AVILA B E F, GAO W, et al. Micro/nanorobots for biomedicine: delivery, surgery, sensing, and detoxification[J]. Science Robotics, 2017, 2(4): eaam6431. DOI: 10.1126/scirobotics.aam6431.
[17] 石振海, 李克智, 李贺军, 等. 航天器热防护材料研究现状与发展趋势[J]. 材料导报, 2017, 21(8): 15-18. (SHI Zhenhai, LI Kezhi, LI Hejun, et al. Research status and application advance of heat resistant materials for space vehicles[J]. Materials Reports, 2007, 21(8): 15-18. (in Chinese))
[18] LI X B, LI L, HU Y J, et al. Bending, buckling and vibration of axially functionally graded beams based on nonlocal strain gradient theory[J]. Composite Structures, 2017, 165: 250-265.
[19] MAHINZARE M, ALIPOUR M J, SADATSAKKAK S A, et al. A nonlocal strain gradient theory for dynamic modeling of a rotary thermo piezo electrically actuated nano FG circular plate[J]. Mechanical Systems and Signal Processing, 2019, 115: 323-337.
[20] KARAMI B, JANGHORBAN M. A new size-dependent shear deformation theory for free vibration analysis of functionally graded/anisotropic nanobeams[J]. Thin-Walled Structures, 2019, 143: 106227.
[21] CHEN P J, PENG J, ZHAO Y C, et al. Prediction of the adhesive behavior of bio-inspired functionally graded materials against rough surfaces[J]. AIP Advances, 2014, 4(6): 067143.
[22] EBRAHIMI F, SHAFIEI N. Application of Eringen’s nonlocal elasticity theory for vibration analysis of rotating functionally graded nanobeams[J]. Smart Structures and Systems, 2016, 17(5): 837-857.
[23] DANESHMEHR A, RAJABPOOR A, POURDAVOOD M. Stability of size dependent functionally graded nanoplate based on nonlocal elasticity and higher order plate theories and different boundary conditions[J]. International Journal of Engineering Science, 2014, 82: 84-100.
[24] BARATI M R, SHAHVERDI H. An analytical solution for thermal vibration of compositionally graded nanoplates with arbitrary boundary conditions based on physical neutral surface position[J]. Mechanics of Advanced Materials and Structures, 2017, 24(10): 840-853.
[25] THAI C H, TRAN T D, PHUNG-VAN P. A size-dependent moving Kriging meshfree model for deformation and free vibration analysis of functionally graded carbon nanotube-reinforced composite nanoplates[J]. Engineering Analysis With Boundary Elements, 2020, 115: 52-63.
[26] ZHANG D G, ZHOU Y H. A theoretical analysis of FGM thin plates based on physical neutral surface[J]. Computational Materials Science, 2008, 44(2): 716-720.