轴向载荷下功能梯度材料Timoshenko梁动力屈曲分析

黄悦; 韩志军; 路国运

doi:10.11858/gywlxb.20180509

轴向载荷下功能梯度材料Timoshenko梁动力屈曲分析

doi: 10.11858/gywlxb.20180509

黄悦^1,,
韩志军^1,*, ,,
路国运²

1.
太原理工大学力学学院, 山西太原 030024
2.
太原理工大学建筑与土木工程学院, 山西太原 030024

基金项目:

国家自然科学基金 11372209

山西省研究生教育改革研究课题 2015JG41

详细信息

作者简介:
黄悦(1991-), 女, 硕士研究生, 主要从事非线性动力屈曲研究.E-mail:huangyue_tyut@163.com

通讯作者:
韩志军(1964-), 男, 博士, 教授, 主要从事非线性动力屈曲研究.E-mail:13073578705@126.com

中图分类号: O344.1;TB333
计量
- 文章访问数: 6980
- HTML全文浏览量: 2936
- PDF下载量: 184
出版历程
- 收稿日期: 2018-01-22
- 修回日期: 2018-02-06

Dynamic Buckling of Functionally Graded Timoshenko Beam under Axial Load

HUANG Yue^1
,,
HAN Zhijun^{1,*
, ,},
LU Guoyun²

1.
College of Mechanics, Taiyuan University of Technology, Taiyuan 030024, China
2.
College of Architecture and Civil Engineering, Taiyuan University of Technology, Taiyuan 030024, China

摘要

摘要: 假设功能梯度材料Timoshenko梁各项物性参数只沿厚度方向按幂函数进行连续变化，研究了功能梯度材料Timoshenko梁的动力屈曲。基于一阶剪切理论，采用Hamilton原理推导出轴向载荷作用下，功能梯度材料Timoshenko梁动力屈曲的控制方程。利用里兹法与棣莫弗公式相结合，获得了功能梯度材料Timoshenko梁在夹支-固支边界条件下动力屈曲临界载荷的解析表达式和屈曲解。应用MATLAB编程计算，讨论了功能梯度材料Timoshenko梁的几何尺寸、梯度指数、模态数、材料构成、泊松比以及弹性模量对临界载荷的影响。结果表明：功能梯度材料Timoshenko梁动力屈曲临界载荷随梁长度的增大而减小，随着梯度指数的增大而减小，随模态数的增大而增大，说明冲击载荷越大，高阶模态越容易被激发；随着泊松比和弹性模量的增大而增大，且泊松比的影响较小，而弹性模量的影响较大。由于剪切项的影响，临界载荷-临界长度的关系曲线在加载端变化趋势平缓。随着模态数的增大，梁的屈曲模态越为复杂。
- 功能梯度材料 /
- 动力屈曲 /
- 临界载荷 /
- 里兹法 /
- 棣莫弗公式
Abstract: In this study, we investigated the dynamic buckling of the functionally graded Timoshenko beam whose property parameters continuously change according to the power function along the thickness direction.Based on the first order shear deformation theory, we derived the governing equation of the dynamic buckling of functionally graded material Timoshenko beams under axial step loading by using the Hamilton's principle.Using the Ritz method combining with the de Moivre's formula, we obtained the buckling solution and the expression of the critical load of the dynamic buckling of functionally graded material Timoshenko beam under the clamped-fixed boundary condition.Then, the influence of geometric size, gradient index, modal number, material composition, Poisson's ratio and elastic modulus on the critical load by MATLAB calculation was discussed.The results show that the critical load of the functionally graded material Timoshenko beam decreases with the increase of beam length and the gradient index, and increases with the increase of the modal number, showing that the higher modal number is more easily excited by the increase of impact load.Furthermore, the critical load increases with the increase of the Poisson's ratio and the elastic modulus, and the effect of elastic modulus is greater than Poisson's ratio.The critical load-critical length curve tends to be gentle at the loading end because of the influence of shear term.Buckling mode of beam becomes more complicated when the modal number increases.
- functionally graded material /
- dynamic buckling /
- critical load /
- Ritz method /
- de Moivre's formula

HTML全文

功能梯度材料是一种集各组分于一体的新型材料，其微观结构和综合性能都呈连续性变化，从而降低了材料的应力集中，提高了材料的性能与功能。因此，功能梯度材料在航天航空、船舶工业、汽车工业、能源工业以及生物医学工业等工程领域有广泛的应用前景，国内外众多学者对功能梯度材料进行了深入的分析和探究^[1-3]。Kocaturk和Akbaş^[4-5]考虑几何非线性，使用有限元法对不同边界条件下功能梯度材料Timoshenko梁在热载作用下的后屈曲进行了分析。Paul^[6]对面内热负荷作用下功能梯度材料Timoshenko梁的后屈曲进行了研究。Şimşek^[7]基于Timoshenko梁理论，首次对功能梯度材料在双向递变情况下的梁动力屈曲进行了讨论。Eltaher等^[8]基于Timoshenko梁理论和最小势能原理，运用有限元法对功能梯度材料Timoshenko梁的动力屈曲进行了分析。Rahimi等^[9]基于Timoshenko梁理论，采用Von-Karman非线性应变-位移关系，对功能梯度材料梁的后屈曲进行了研究。Li^[10]比较了轴向载荷作用下功能梯度材料Timoshenko梁和欧拉-伯努利梁的临界屈曲载荷之间的关系。Kahya等^[11]采用有限元软件对功能梯度梁的动力屈曲以及自由振动进行了分析。Kiani^[12]基于欧拉-伯努利梁理论，对不同热载荷作用下功能梯度材料梁的动力屈曲进行了研究。Anandrao^[13]运用里兹法和有限元法对均匀细长功能梯度梁热后屈曲进行了讨论。Akbaş等^[14]利用Newton-Raphson迭代法结合有限元法研究了热载荷作用下功能梯度梁的后屈曲行为。Rychlewska^[15]通过等效梁近似功能梯度梁，对功能梯度梁进行了动力屈曲分析。综上所述，功能梯度梁的屈曲研究主要以数值计算为主。

本研究基于一阶剪切理论，利用Hamilton原理推导出功能梯度材料Timoshenko梁动力屈曲的控制方程；采用里兹法结合棣莫弗公式，导出功能梯度材料Timoshenko梁在夹支-固支边界条件下动力屈曲临界载荷的解析表达式；应用MATLAB软件编程计算，讨论功能梯度材料Timoshenko梁的几何尺寸、梯度指数、模态数、材料构成、泊松比以及弹性模量对临界载荷的影响。

1. 功能梯度材料物性参数

功能梯度材料Timoshenko梁的各项物性参数只沿厚度(z)方向按幂函数连续变化，表示为^[16]

$弹性模量~~~~{E_z} = \left( {{E_{\rm{c}}} - {E_{\rm{m}}}} \right){\left( {\frac{1}{2} + \frac{z}{h}} \right)^k} + {E_{\rm{m}}}$

(1)

$密度~~~~{\rho _z} = \left( {{\rho _{\rm{c}}} - {\rho _{\rm{m}}}} \right){\left( {\frac{1}{2} + \frac{z}{h}} \right)^k} + {\rho _{\rm{m}}}$

(2)

$泊松比~~~~{\mu _z} = \left( {{\mu _{\rm{c}}} - {\mu _{\rm{m}}}} \right){\left( {\frac{1}{2} + \frac{z}{h}} \right)^k} + {\mu _{\rm{m}}}$

(3)

$剪切模量~~~~{G_z} = \frac{{{E_z}}}{{2\left( {1 - {\mu _z}} \right)}}$

(4)

式中：E_c是陶瓷的弹性模量，E_m是金属的弹性模量，ρ_c是陶瓷的密度，ρ_m是金属的密度，μ_c是陶瓷的泊松比，μ_m是金属的泊松比, k是梯度指数。

2. 动力控制方程的推导

如图 1所示，功能梯度材料Timoshenko梁的长为l，宽为b，高为h，材料各项性质只沿高度方向呈连续性变化，剖面中心线为x轴，设其轴向位移为u(x, t)，横向位移为w(x, t)，横截面转角为φ。基于平截面假定，由一阶剪切理论，Timoshenko梁的位移及应变^[17]可表示为

$\left\{ \begin{array}{l} u = {u_0} - z\varphi \\ w = {w_0} \end{array} \right.$

(5)

$\left\{ \begin{array}{l} {\varepsilon _x} = \varepsilon _x^0 + z{k_x}\\ {\gamma _{xz}} = \gamma _{xz}^0 \end{array} \right.$

(6)

$其中~~~~\varepsilon _x^0 = \frac{{\partial {u_0}}}{{\partial x}},\;\;\;\;\gamma _{xz}^0 = \frac{{\partial {w_0}}}{{\partial x}} - \varphi ,\;\;\;\;{k_x} = - \frac{{\partial \varphi }}{{\partial x}}$

(7)

图 1 受轴向载荷作用的Timoshenko梁

Figure 1. Timoshenko beam under axial load

下载: 全尺寸图片幻灯片

式中：u₀是杆的中性面在x轴方向的位移，w₀是杆的中性面在z轴方向的位移，ε_x⁰是杆的中性面在x方向的正应变，γ_xz⁰是杆的中性面在xz方向的切应变，k_x是杆件任意点在x轴方向上的曲率。

其内力、剪力及内力矩可表示为

$内力~~~~{N_x} = {A_{11}}\varepsilon _x^0 + {B_{11}}{k_x}$

(8)

$剪力~~~~{Q_x} = K{A_{55}}\left( {\frac{{\partial {w_0}}}{{\partial x}} - \varphi } \right)$

(9)

$内力矩~~~~{M_x} = {B_{11}}\varepsilon _x^0 + {D_{11}}{k_x}$

(10)

式中：A₁₁为拉伸刚度系数; B₁₁为耦合刚度系数; D₁₁为扭转刚度系数; A₅₅为剪切刚度系数; K为剪切修正系数，取K=5/6。

${A_{11}} = \int_{ - h/2}^{h/2} {\frac{{{E_z}}}{{1 - \mu _z^2}}{\rm{d}}z}$

(11)

${B_{11}} = \int_{ - h/2}^{h/2} {\frac{{z{E_z}}}{{1 - \mu _z^2}}{\rm{d}}z}$

(12)

${D_{11}} = \int_{ - h/2}^{h/2} {\frac{{{z^2}{E_z}}}{{1 - \mu _z^2}}{\rm{d}}z}$

(13)

${A_{55}} = \int_{ - h/2}^{h/2} {\frac{{{E_z}}}{{2\left( {1 + {\mu _z}} \right)}}{\rm{d}}z}$

(14)

$应变能~~~~U = \frac{1}{2}\int_0^l {\int_0^b {\int_{ - h/2}^{h/2} {\left( {{N_x}\varepsilon _x^0 + {M_x}{k_x} + {Q_x}{\gamma _{xz}}} \right){\rm{d}}x{\rm{d}}y{\rm{d}}z} } }$

(15)

$动能~~~~T = \frac{1}{2}\int_0^l {\int_0^b {\int_{ - h/2}^{h/2} {{\rho _z}\left[ {{{\left( {\frac{{\partial u}}{{\partial t}}} \right)}^2} + {{\left( {\frac{{\partial w}}{{\partial t}}} \right)}^2}} \right]{\rm{d}}x{\rm{d}}y{\rm{d}}z} } }$

(16)

$外力功~~~~W = \frac{1}{2}\int_0^l {\int_0^b {\int_{ - h/2}^{h/2} {{N_x}{{\left( {\frac{{\partial w}}{{\partial x}}} \right)}^2}{\rm{d}}x{\rm{d}}y{\rm{d}}z} } }$

(17)

定义广义变量为

$\left( {{I_0},{I_1},{I_2}} \right) = \int_{ - h/2}^{h/2} {{\rho _z}\left( {1,z,{z^2}} \right){\rm{d}}z}$

(18)

将(15)式~(18)式代入Hamilton原理，即 $\delta \int_{{t_1}}^{{t_2}} {\left( {T{\rm{ - }}U + W} \right){\rm{d}}t = 0}$ ，由u₀、v₀、φ的变分系数取零可得

${\rm{ \mathsf{ δ} }}{u_0}: - {I_0}\frac{{{\partial ^2}{u_0}}}{{\partial {t^2}}} + {I_1}\frac{{{\partial ^2}\varphi }}{{\partial {t^2}}} + {A_{11}}\frac{{{\partial ^2}{u_0}}}{{\partial {x^2}}} - {B_{11}}\frac{{{\partial ^2}\varphi }}{{\partial {x^2}}} = 0$

(19)

${\rm{ \mathsf{ δ} }}{w_0}: - {I_0}\frac{{{\partial ^2}{w_0}}}{{\partial {t^2}}} - N\left( t \right)\frac{{{\partial ^2}{w_0}}}{{\partial {x^2}}} + K{A_{55}}\frac{{{\partial ^2}{w_0}}}{{\partial {x^2}}} - K{A_{55}}\frac{{\partial \varphi }}{{\partial x}} = 0$

(20)

${\rm{ \mathsf{ δ} }}\varphi :{I_1}\frac{{{\partial ^2}{u_0}}}{{\partial {t^2}}} + {I_2}\frac{{{\partial ^2}\varphi }}{{\partial {t^2}}} - {B_{11}}\frac{{{\partial ^2}{u_0}}}{{\partial {x^2}}} + {D_{11}}\frac{{{\partial ^2}\varphi }}{{\partial {x^2}}} - K{A_{55}}\varphi + K{A_{55}}\frac{{\partial {w_0}}}{{\partial x}} = 0$

(21)

不考虑轴向惯性^[18]，将(21)式对x求导可得

$- {I_2}\frac{{{\partial ^3}\varphi }}{{\partial {t^2}\partial x}} + {D_{11}}\frac{{{\partial ^3}\varphi }}{{\partial {x^3}}} - K{A_{55}}\frac{{\partial \varphi }}{{\partial x}} + K{A_{55}}\frac{{{\partial ^2}{w_0}}}{{\partial {x^2}}} = 0$

(22)

对(20)式进行偏导计算，得

$\left\{ \begin{array}{l} K{A_{55}}\frac{{\partial \varphi }}{{\partial x}} = - {I_0}\frac{{{\partial ^2}{w_0}}}{{\partial {t^2}}} - N\left( t \right)\frac{{{\partial ^2}{w_0}}}{{\partial {x^2}}} + K{A_{55}}\frac{{{\partial ^2}{w_0}}}{{\partial {x^2}}}\\ K{A_{55}}\frac{{{\partial ^3}\varphi }}{{\partial x\partial {t^2}}} = - {I_0}\frac{{{\partial ^4}{w_0}}}{{\partial {t^4}}} - N\left( t \right)\frac{{{\partial ^4}{w_0}}}{{\partial {x^2}\partial {t^2}}} + K{A_{55}}\frac{{{\partial ^4}{w_0}}}{{\partial {x^2}\partial {t^2}}}\\ K{A_{55}}\frac{{{\partial ^2}\varphi }}{{\partial {x^2}}} = - {I_0}\frac{{{\partial ^3}{w_0}}}{{\partial {t^2}\partial x}} - N\left( t \right)\frac{{{\partial ^3}{w_0}}}{{\partial {x^3}}} + K{A_{55}}\frac{{{\partial ^3}{w_0}}}{{\partial {x^3}}}\\ K{A_{55}}\frac{{{\partial ^3}\varphi }}{{\partial {x^3}}} = - {I_0}\frac{{{\partial ^4}{w_0}}}{{\partial {t^2}\partial {x^2}}} - N\left( t \right)\frac{{{\partial ^4}{w_0}}}{{\partial {x^4}}} + K{A_{55}}\frac{{{\partial ^4}{w_0}}}{{\partial {x^4}}} \end{array} \right.$

(23)

将(23)式代入(22)式，得到控制方程为

$\begin{array}{l} \frac{{{D_{11}}}}{{K{A_{55}}}}\left[ { - {I_0}\frac{{{\partial ^4}w}}{{\partial {x^2}\partial {t^2}}} - N\left( t \right)\frac{{{\partial ^4}w}}{{\partial {x^4}}}} \right] + {D_{11}}\frac{{{\partial ^4}w}}{{\partial {x^4}}} + \frac{{{I_2}}}{{K{A_{55}}}}\left[ {{I_0}\frac{{{\partial ^4}w}}{{\partial {t^4}}} + N\left( t \right)\frac{{{\partial ^4}w}}{{\partial {x^2}\partial {t^2}}}} \right] - \\ {I_2}\frac{{{\partial ^4}w}}{{\partial {x^2}\partial {t^2}}} + {I_0}\frac{{{\partial ^2}w}}{{\partial {t^2}}} + N\left( t \right)\frac{{{\partial ^2}w}}{{\partial {x^2}}} = 0 \end{array}$

(24)

3. 动力屈曲临界载荷

3.1 用里兹法化简控制方程

设一端夹支、一端固支梁的动力屈曲解的形式为^[19]

$w = T\left( t \right)\left( {A\sin \frac{{{n_1}{\rm{ \mathsf{ π} }}x}}{l} + B\sin \frac{{{n_2}{\rm{ \mathsf{ π} }}x}}{l}} \right)$

(25)

式中：n₁、n₂为屈曲的模态阶数，A、B为待求系数。动力屈曲解要满足下列一端夹支、一端固支梁的边界条件

$w\left( 0 \right) = \frac{{\partial w}}{{\partial x}}\left| {_{x = 0}} \right. = w\left( l \right) = \frac{{\partial w}}{{\partial x}}\left| {_{x = l}} \right. = 0$

(26)

将(26)式代入(25)式可得

$\left\{ \begin{array}{l} A\frac{{{n_1}{\rm{ \mathsf{ π} }}}}{l} + B\frac{{{n_2}{\rm{ \mathsf{ π} }}}}{l} = 0\\ A\frac{{{n_1}{\rm{ \mathsf{ π} }}}}{l}\cos \left( {{n_1}{\rm{ \mathsf{ π} }}} \right) + B\frac{{{n_2}{\rm{ \mathsf{ π} }}}}{l}\cos \left( {{n_2}{\rm{ \mathsf{ π} }}} \right) = 0 \end{array} \right.$

(27)

若(27)式有解则cos(n₁π)与cos(n₂π)奇偶性相同，即n₂=n₁+2。又 $A\frac{{{n_1}{\rm{ \mathsf{ π} }}}}{l} + B\frac{{{n_2}{\rm{ \mathsf{ π} }}}}{l} = 0$ ，即 $\frac{A}{B} = - \frac{{{n_2}}}{{{n_1}}}$ 。取 $A = 1, B = - \frac{{{n_1}}}{{{n_2}}}$ 。则梁的动力屈曲解的形式化为

$w = T\left( t \right)\left( {\sin \frac{{{n_1}{\rm{ \mathsf{ π} }}x}}{l} - \frac{{{n_1}}}{{{n_2}}}\sin \frac{{{n_2}{\rm{ \mathsf{ π} }}x}}{l}} \right)$

(28)

式中：n₁、n₂是屈曲的模态阶数，且n₁=1, 2, 3，…; n₂=n₁+2。

对(28)式求偏导可得

$\left\{ \begin{array}{l} {w_{xx}} = T\frac{{{n_1}{{\rm{ \mathsf{ π} }}^2}}}{{{l^2}}}\left( { - {n_1}\sin \frac{{{n_1}{\rm{ \mathsf{ π} }}x}}{l} + {n_2}\sin \frac{{{n_2}{\rm{ \mathsf{ π} }}x}}{l}} \right)\\ {w_{tt}} = {T_{tt}}\left( {\sin \frac{{{n_1}{\rm{ \mathsf{ π} }}x}}{l} - \frac{{{n_1}}}{{{n_2}}}\sin \frac{{{n_2}{\rm{ \mathsf{ π} }}x}}{l}} \right)\\ {w_{xxtt}} = {T_{tt}}\frac{{{n_1}{{\rm{ \mathsf{ π} }}^2}}}{{{l^2}}}\left( { - {n_1}\sin \frac{{{n_1}{\rm{ \mathsf{ π} }}x}}{l} + {n_2}\sin \frac{{{n_2}{\rm{ \mathsf{ π} }}x}}{l}} \right)\\ {w_{xxxx}} = T\frac{{{n_1}{{\rm{ \mathsf{ π} }}^4}}}{{{l^4}}}\left( {n_1^3\sin \frac{{{n_1}{\rm{ \mathsf{ π} }}x}}{l} - n_2^3\sin \frac{{{n_2}{\rm{ \mathsf{ π} }}x}}{l}} \right)\\ {w_{tttt}} = {T_{tttt}}\left( {\sin \frac{{{n_1}{\rm{ \mathsf{ π} }}x}}{l} - \frac{{{n_1}}}{{{n_2}}}\sin \frac{{{n_2}{\rm{ \mathsf{ π} }}x}}{l}} \right) \end{array} \right.$

(29)

3.2 由棣莫弗公式求解临界载荷

由棣莫弗公式 ${\left( {{\rm{cos}}\theta + {\rm{isin}}\theta } \right)^n} = {\rm{cos}}\left( {n\theta } \right) + {\rm{isin}}\left( {n\theta } \right)$ 可得

$\frac{{\sin \frac{{n{\rm{ \mathsf{ π} }}x}}{l}}}{{\sin \frac{{{\rm{ \mathsf{ π} }}x}}{l}}} = C_n^1{\left( {\cos \frac{{{\rm{ \mathsf{ π} }}x}}{l}} \right)^{n - 1}} - C_n^3{\left( {\cos \frac{{{\rm{ \mathsf{ π} }}x}}{l}} \right)^{n - 3}}{\left( {\sin \frac{{{\rm{ \mathsf{ π} }}x}}{l}} \right)^2} + \cdots$

(30)

将(29)式除以 ${\rm{sin}}\frac{{{\rm{ \mathsf{ π} }}x}}{l}$ ，当梁发生屈曲，即x=l时有

$\left\{ \begin{array}{l} \frac{{{w_{xx}}}}{{\sin \frac{{{\rm{ \mathsf{ π} }}x}}{l}}} = T\frac{{{n_1}{{\rm{ \mathsf{ π} }}^2}}}{{{l^2}}}\left( { - n_1^2 + n_2^2} \right)\\ \frac{{{w_{tt}}}}{{\sin \frac{{{\rm{ \mathsf{ π} }}x}}{l}}} = {T_{tt}}\left( {{n_1} - {n_1}} \right) = 0\\ \frac{{{w_{xxtt}}}}{{\sin \frac{{{\rm{ \mathsf{ π} }}x}}{l}}} = {T_{tt}}\frac{{{n_1}{{\rm{ \mathsf{ π} }}^2}}}{{{l^2}}}\left( { - n_1^2 + n_2^2} \right)\\ \frac{{{w_{xxxx}}}}{{\sin \frac{{{\rm{ \mathsf{ π} }}x}}{l}}} = T\frac{{{n_1}{{\rm{ \mathsf{ π} }}^4}}}{{{l^4}}}\left( {n_1^4 - n_2^4} \right)\\ \frac{{{w_{tttt}}}}{{\sin \frac{{{\rm{ \mathsf{ π} }}x}}{l}}} = {T_{tttt}}\left( {{n_1} - {n_1}} \right) = 0 \end{array} \right.$

(31)

将(31)式代入控制方程得

${T_{tt}} + T\left( {\alpha - \beta } \right) = 0$

(32)

其中

$\alpha = N\left( t \right)\frac{{{D_{11}}\left( {n_1^2 + n_2^2} \right){{\rm{ \mathsf{ π} }}^2} + K{A_{55}}{l^2}}}{{\left[ { - {D_{11}}{I_0} + {I_2}N\left( t \right) - {I_2}K{A_{55}}} \right]}},\;\;\;\;\beta = \frac{{{D_{11}}\left( {n_1^2 + n_2^2} \right){{\rm{ \mathsf{ π} }}^2}K{A_{55}}}}{{\left[ { - {D_{11}}{I_0} + {I_2}N\left( t \right) - {I_2}K{A_{55}}} \right]{l^2}}}$

当α-β＜0时, 方程的解为指数形式的发散解，表明系统经受微小的扰动后，运动为无界，因此系统呈现不稳定的状态。当α-β＞0时，表明系统围绕平衡位置做微小振动，因此系统呈现稳定平衡的状态^[19]。当α-β=0时为Timoshenko梁屈曲临界状态，由此推出动力屈曲临界载荷的表达式为

${N_{{\rm{cr}}}} = \frac{{{D_{11}}\left( {n_1^2 + n_2^2} \right){{\rm{ \mathsf{ π} }}^2}K{A_{55}}}}{{{D_{11}}\left( {n_1^2 + n_2^2} \right){{\rm{ \mathsf{ π} }}^2} + K{A_{55}}l_{{\rm{cr}}}^2}}$

(33)

(33)式表明：剪切项对动力屈曲临界载荷有一定的影响，当临界长度l_cr=0时，临界载荷为KA₅₅。

4. 数值计算

运用MATLAB软件对(33)式进行编程计算，材料选取陶瓷/钛功能梯度材料、陶瓷/铁功能梯度材料和陶瓷/铜功能梯度材料。梁的宽度b=2cm，厚度h=3cm。讨论梁厚度、模态数(n)、梯度指数(k)、泊松比(μ)、弹性模量(E)以及材料构成对临界载荷的影响。材料各参数如表 1所示。

表 1 材料各项参数

Table 1. Material parameters

Material	Elasticity modulus E/GPa	Density ρ/ (g·cm^-3)	Poisson's ratio μ
Ceramic	385	3.96	0.23
Titanium	108.5	4.54	0.41
Iron	155	7.86	0.291
Copper	119	8.96	0.326

下载: 导出CSV

| 显示表格

图 2所示为n=1、k=1/2时，不同材料构成下临界载荷N_cr与临界长度l_cr的关系曲线。当h、n、k、μ以及E一定，材料构成不同时，其N_cr随着l_cr的增加而减小。当l_cr相同时，陶瓷/铁材料的临界载荷值高于陶瓷-钛材料，且陶瓷/钛材料梁与陶瓷/铜材料梁的N_cr-l_cr曲线在l_cr超过1m时接近重合。由于剪切项的影响，曲线在加载端变化趋势平缓。

图 2 不同材料构成下临界载荷与临界长度的关系

Figure 2. Relationship between critical load and critical length for different materials

下载: 全尺寸图片幻灯片

图 3所示为陶瓷/钛功能梯度材料梁在n=1时，不同k下N_cr与l_cr的关系曲线。当梁的h、n、μ、E以及材料构成一定，k不同时，其N_cr随着l_cr的增大而减小。当l_cr相同时，N_cr随着的k的增大而减小。由于剪切项的影响，曲线在加载端变化趋势也平缓，如图 2所示。

图 3 不同梯度指数下临界载荷与临界长度的关系

Figure 3. Relationship between critical load and critical length for different gradient indexes

下载: 全尺寸图片幻灯片

图 4所示为陶瓷/钛功能梯度材料梁在k=1/2时，不同n下N_cr与l_cr的关系曲线。当梁的h、k、μ、E以及材料构成一定，n不同时，其N_cr随着l_cr的增大而减小。当l_cr相同时，N_cr随着n的增大而增大。这表明冲击载荷越大，高阶模态越容易被激发。

图 4 不同模态数下临界载荷与临界长度的关系

Figure 4. Relationship between critical load and critical length for different modal numbers

下载: 全尺寸图片幻灯片

图 5表示陶瓷/钛功能梯度材料梁在k=1/2、n=1时，不同h下N_cr与l_cr的关系曲线。当梁的n、k、μ、E以及材料构成一定，h不同时，其N_cr随着l_cr的增大而减小。当l_cr相同时，N_cr随着h的增大而增大。

图 5 不同截面高度下临界载荷与临界长度的关系

Figure 5. Relationship between critical load and critical length for different heights

下载: 全尺寸图片幻灯片

图 6表示k=1/2、n=1、μ=0.41时，不同E下N_cr与l_cr的关系曲线。当梁的h、k、μ、n一定，E不同时，其N_cr随着l_cr的增大而减小。当l_cr相同时，N_cr随着E的增大而增大。

图 6 不同弹性模量下临界载荷与临界长度的关系

Figure 6. Relationship between critical load and critical length for different elasticity modulus

下载: 全尺寸图片幻灯片

图 7表示k=1/3、n=1、E=100GPa时，不同μ下N_cr与l_cr的关系曲线。当梁的h、k、E、n一定，μ不同时，其N_cr随着l_cr的增大而减小。当l_cr相同时，不同μ下曲线接近重合，表明μ对N_cr的影响不大。

图 7 不同泊松比下临界载荷与临界长度的关系

Figure 7. Relationship between critical load and critical length for different Poisson's ratios

下载: 全尺寸图片幻灯片

图 8~图 11是陶瓷/钛功能梯度材料梁在k=1/2时，不同模态数下的屈曲模态图。可见，随着模态数n的增大，梁的屈曲模态越发复杂。固定端屈曲模态的幅值较大，对比不同模态下的固支端与夹支端，其幅值比例范围为1.62~2.88。

图 8 梁屈曲模态数取n=1模态图

Figure 8. Buckling mode of beam for n=1

下载: 全尺寸图片幻灯片

图 9 梁屈曲模态数取n=2模态图

Figure 9. Buckling mode of beam for n=2

下载: 全尺寸图片幻灯片

图 10 梁屈曲模态数取n=3模态图

Figure 10. Buckling mode of beam for n=3

下载: 全尺寸图片幻灯片

图 11 梁屈曲模态数取n=4模态图

Figure 11. Buckling mode of beam for n=4

下载: 全尺寸图片幻灯片

5. 结论

(1) 基于一阶剪切理论，通过Hamilton原理推导出功能梯度材料Timoshenko梁的动力屈曲控制方程。

(2) 用里兹法与棣莫弗公式相结合，得到了功能梯度材料Timoshenko梁在夹支-固支边界条件下动力屈曲临界载荷的解析表达式及屈曲解。

(3) 应用MATLAB软件，讨论了功能梯度材料Timoshenko梁的几何尺寸、梯度指数、模态数、材料构成、泊松比以及弹性模量对临界载荷的影响。结果表明：功能梯度材料Timoshenko梁的动力屈曲临界载荷随着梁长度、梯度指数的增大而减小；临界载荷随模态数的增大而增大，且冲击载荷越大，高阶模态越容易被激发；临界载荷随泊松比及弹性模量的增大而增大，且泊松比的影响不大，而弹性模量的影响较大。由于剪切项的影响，临界载荷-临界长度的关系曲线在加载端变化趋势平缓；随着模态数的增大，梁的屈曲模态越发复杂。

图 1 受轴向载荷作用的Timoshenko梁

Figure 1. Timoshenko beam under axial load

下载: 全尺寸图片幻灯片

图 2 不同材料构成下临界载荷与临界长度的关系

Figure 2. Relationship between critical load and critical length for different materials

下载: 全尺寸图片幻灯片

图 3 不同梯度指数下临界载荷与临界长度的关系

Figure 3. Relationship between critical load and critical length for different gradient indexes

下载: 全尺寸图片幻灯片

图 4 不同模态数下临界载荷与临界长度的关系

Figure 4. Relationship between critical load and critical length for different modal numbers

下载: 全尺寸图片幻灯片

图 5 不同截面高度下临界载荷与临界长度的关系

Figure 5. Relationship between critical load and critical length for different heights

下载: 全尺寸图片幻灯片

图 6 不同弹性模量下临界载荷与临界长度的关系

Figure 6. Relationship between critical load and critical length for different elasticity modulus

下载: 全尺寸图片幻灯片

图 7 不同泊松比下临界载荷与临界长度的关系

Figure 7. Relationship between critical load and critical length for different Poisson's ratios

下载: 全尺寸图片幻灯片

图 8 梁屈曲模态数取n=1模态图

Figure 8. Buckling mode of beam for n=1

下载: 全尺寸图片幻灯片

图 9 梁屈曲模态数取n=2模态图

Figure 9. Buckling mode of beam for n=2

下载: 全尺寸图片幻灯片

图 10 梁屈曲模态数取n=3模态图

Figure 10. Buckling mode of beam for n=3

下载: 全尺寸图片幻灯片

图 11 梁屈曲模态数取n=4模态图

Figure 11. Buckling mode of beam for n=4

下载: 全尺寸图片幻灯片

表 1 材料各项参数

Table 1. Material parameters

Material	Elasticity modulus E/GPa	Density ρ/ (g·cm^-3)	Poisson's ratio μ
Ceramic	385	3.96	0.23
Titanium	108.5	4.54	0.41
Iron	155	7.86	0.291
Copper	119	8.96	0.326

下载: 导出CSV

参考文献(19)

[1]	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
[2]	GUPTA A, TALHA M.Recent development in modeling and analysis of functionally graded materials and structures[J]. Progress in Aerospace Sciences, 2015, 79:1-14. doi: 10.1016/j.paerosci.2015.07.001
[3]	NAEBE M, SHIRVANIMOGHADDAM K.Functionally graded materials:a review of fabrication and properties[J]. Applied Materials Today, 2016, 5:223-245. doi: 10.1016/j.apmt.2016.10.001
[4]	KOCATURK T, AKBAŞ Ş D.Post-buckling analysis of Timoshenko beams made of functionally graded material under thermal loading[J]. Structural Engineering & Mechanics, 2012, 41(6):775-789. doi: 10.1007/s10999-010-9132-4
[5]	KOCATURK T, AKBAŞ Ş D.Post-buckling analysis of Timoshenko beams with various boundary conditions under non-uniform thermal loading[J]. Structural Engineering & Mechanics, 2011, 40(3):347-371. http://www.koreascience.or.kr/article/ArticleFullRecord.jsp?cn=KJKHB9_2011_v40n3_347
[6]	PAUL A.Non-linear thermal post-buckling analysis of FGM Timoshenko beam under non-uniform temperature rise across thickness[J]. Engineering Science & Technology:An International Journal, 2016, 19(3):1608-1625. https://www.sciencedirect.com/science/article/pii/S2215098616302853
[7]	ŞIMŞEK M.Buckling of Timoshenko beams composed of two-dimensional functionally graded material (2D-FGM) having different boundary conditions[J]. Composite Structures, 2016, 149:304-314. doi: 10.1016/j.compstruct.2016.04.034
[8]	ELTAHER M A, KHAIRY A, SADOUN A M, et al.Static and buckling analysis of functionally graded Timoshenko nanobeams[J]. Applied Mathematics &Computation, 2014, 229(229):283-295. https://www.sciencedirect.com/science/article/pii/S0096300313013556
[9]	RAHIMI G H, GAZOR M S, HEMMATNEZHAD M, et al.On the postbuckling and free vibrations of FG Timoshenko beams[J]. Composite Structures, 2013, 95(1):247-253. https://www.sciencedirect.com/science/article/pii/S0263822312003625
[10]	LI S R.Relations between buckling loads of functionally graded Timoshenko and homogeneous Euler-Bernoulli beams[J]. Composite Structures, 2013, 95(1):5-9.
[11]	KAHYA V, TURAN M.Finite element model for vibration and buckling of functionally graded beams based on the first-order shear deformation theory[J]. Composites Part B:Engineering, 2017, 109:108-115. doi: 10.1016/j.compositesb.2016.10.039
[12]	KIANI Y.Thermal buckling analysis of functionally graded material beams[J]. International Journal of Mechanics & Materials in Design, 2010, 6(3):229-238. doi: 10.1007/s10999-010-9132-4
[13]	ANANDRAO K S.Thermal post-buckling analysis of uniform slender functionally graded material beams[J]. Structural Engineering & Mechanics, 2010, 36(5):545-560. https://www.researchgate.net/publication/264077187_Thermal_post-buckling_analysis_of_uniform_slender_functionally_graded_material_beams
[14]	AKBAŞ Ş D, KOCATÜRK T.Post-buckling analysis of functionally graded three-dimensional beams under the influence of temperature[J]. Journal of Thermal Stresses, 2013, 36(12):1233-1254. doi: 10.1080/01495739.2013.788397
[15]	RYCHLEWSKA J.A new approach for buckling analysis of axially functionally graded beams[J]. Journal of Applied Mathematics & Computational Mechanics, 2015, 14(2):95-102.
[16]	FEREZQI H Z, TAHANI M, TOUSSI H E.Analytical approach to free vibrations of cracked Timoshenko beams made of functionally graded materials[J]. Mechanics of Advanced Materials & Structures, 2010, 17(5):353-365. doi: 10.1080/15376494.2010.488608?scroll=top&needAccess=true
[17]	徐华, 李世荣.一阶剪切理论下功能梯度梁与均匀梁静态解之间的相似关系[J].工程力学, 2012, 29(4):161-167. http://www.cnki.com.cn/Article/CJFDTotal-GTLX201403013.htm XU H, LI S R.Analogous relationship between the static solutions of functionally graded beams and homogeneous beams based on the first-order shear deformation theory[J]. Engineering Mechanics, 2012, 29(4):161-167. http://www.cnki.com.cn/Article/CJFDTotal-GTLX201403013.htm
[18]	KIRCHHOFF G.Ueber das Gleichgewicht und die Bewegung eines unendlich dünnen elastischen Stabes[J]. Journal FüR Die Reine Und Angewandte Mathematik, 2009, 1859(56):285-313. https://eudml.org/doc/147766
[19]	李楠, 韩志军, 路国运.基于里兹法研究复合材料层合板的动力屈曲问题[J].振动与冲击, 2016, 35(10):180-184. http://industry.wanfangdata.com.cn/jt/Detail/Periodical?id=Periodical_zdycj201610029 LI N, HAN Z J, LU G Y.Research on dynamic buckling of laminated composite plates using Ritz method[J]. Journal of Vibration and Shock, 2016, 35(10):180-184. http://industry.wanfangdata.com.cn/jt/Detail/Periodical?id=Periodical_zdycj201610029