基于改进AMDM的弹性边界变截面Timoshenko梁自由振动分析

Free vibration analysis of non-uniform Timoshenko beams under elastic boundary conditions using improved AMDM

摘要

为改善Adomian修正分解法（adomian modified decomposition method， AMDM）在求解弹性边界变截面Timoshenko梁振动微分方程时的收敛性， 把自适应分段技术引入AMDM的计算流程， 通过使截面形状相关项的Taylor级数积分绝对误差估计符合特定阈值， 自动确定具有不同形状函数幂次和锥度的变截面梁分段长度与数量 . 分析了锥度、跨高比、形状函数幂次以及边界支承刚度对线性、幂函数和指数函数变截面梁固有频率的影响， 并探讨了影响收敛性的因素. 数值仿真显示， 与传统AMDM相比， 计算耗时在不同情形下缩短了46%~91%. 计算结果与有限元仿真相比， 误差不超过0.16%， 验证了本文方法的有效性和准确性. 研究表明， 固有频率对跨高比敏感度从大到小依次为幂函数、线性和指数函数变截面梁；边界支承刚度的强弱决定了固有频率-锥度曲线的增减趋势.

Abstract

To improve the convergence performance of the AMDM （adomian modified decomposition method） in solving the differential equations of vibration problem of non-uniform Timoshenko beams under elastic boundary conditions， the adaptive segmentation technology is introduced to the calculation process of the AMDM. By making the integral absolute error estimation of Taylor expansion of section shape dependent terms meet specific thresholds， the length of each segment and the number of segments are obtained automatically for non-uniform beams with shape functions of different power indices and taper ratios. The effects of taper ratio， length-to-height ratio， power index of shape functions and boundary stiffness on natural frequencies of linear， power， and exponential non-uniform beams are analyzed， and the factors affecting convergence are discussed. The numerical results show that the reduction of calculation time ranges from 46% to 91% compared with conventional AMDM in different cases. In addition， compared with the finite element simulation， the error of the calculation results in this paper is within 0.16%， which verifies the effectiveness and accuracy of the method. It is found that the natural frequency exhibits varying sensitivity to length-to-height ratio successively from large to small for power， linear， and exponential non-uniform beams. The boundary support stiffness determines the increasing or decreasing trend of the frequency-taper curve.

Graphical abstract

关键词

变截面 / Timoshenko梁 / 边界条件 / 改进的AMDM / 形状函数 / 固有频率

Key words

non-uniform sections / Timoshenko beams / boundary conditions / improved AMDM / shape functions / natural frequency

变截面梁广泛地应用于机械、建筑等领域，其自由振动问题受到相关学者的普遍关注.跨高比大的梁的振动分析可以用经典的Euler-Bernoulli梁理论建模^［1-5］.

对于跨高比小的梁而言， 横截面内剪切变形的影响不容忽略. 因此， 必须采用更为精确的Timoshenko梁理论对经典理论进行修正^［6］. 潘旦光等^［7］利用模态摄动法求解得到线性变截面Timoshenko梁的特征值和振型. 葛仁余等^［8］采用插值矩阵法求解了轴向功能梯度线性变截面Timoshenko梁的固有频率. 黄意新等^［9］利用Chebyshev谱方法进行了复杂边界条件下轴向功能梯度线性变截面Timoshenko梁的动力学分析. Molina-Villegas等^［10］给出了轴向非均匀Timoshenko梁振动问题的封闭解.除了上述提到的几种解法， Timoshenko梁自由振动问题的求解方法还包括微分求积单元法^［11］，Lagrangian法^［12］和Rayleigh-Ritz法^［13］.

在工程实践中，受弹性约束的情形比受经典约束的情形更具一般性.Goel^［14］分析了扭转约束对线性锥形梁的作用.Grossi等^［15-16］利用Rayleigh-Ritz法和Rayleigh-Schmidt法研究了弹性约束下锥形梁的振动问题. Ho等^［17］基于微分变换法分析了常规弹性约束下非均匀梁的振动问题. 文献［18］和文献［19］利用Bessel函数法分析了柔性边界条件下锥形梁的自由振动问题. Galerkin法、传递矩阵法、谱方法和Rayleigh-Ritz法要么在构造近似解时常常需要对形函数或偏微分方程积分做近似处理， 要么存在精度和计算效率等问题.

AMDM是美国数学物理学家George Adomian在20世纪80年代提出的求解线性和非线性边值问题近似解的一种通用数学方法， 它将未知函数分解为级数形式， 将积分逆运算算子应用于方程或方程组， 适用于求解线性/非线性微分、时滞/积分微分以及偏微分方程和方程组^［20-21］. 与上述提到的各种方法相比， AMDM是一种具有普适性的高效、高精度计算方法. 目前， 在基于AMDM的振动分析中涉及的结构， 包括变截面Euler-Bernoulli梁^［22-24］和均匀Timoshenko梁^［25］. 郭银凤^［26］指出，对AMDM的改进主要有扩大其应用范围、简化Adomian多项式的计算、提高收敛速度等方向.

AMDM的核心是利用Taylor级数展开逼近原函数. 尽管在较小的区域内可以实现快速收敛， 但在较大区域的收敛速度却相当缓慢，而且截断级数解在较大区域内是不精确的^［27-28］.因此，采用传统的AMDM求解变截面梁的自由振动问题时，迭代运算时间长、计算效率低，特别是针对锥度较大的线性变截面梁以及非线性变截面梁时，常常无法收敛， 使求解过程变得异常困难.

自适应网格划分是有限元仿真获得高精度解答的一种高性能计算方案^［29-30］.本文将自适应分段技术引入AMDM计算过程，采用积分绝对误差（integral absolute error， IAE）作为误差指标^［31］，通过逐渐缩小分段长度让截面形状相关项的Taylor展开式与原函数的IAE满足特定阈值， 由此确定每段的长度和分段的数量. 采用改进的AMDM求解各种弹性边界条件下变截面Timoshenko梁自由振动问题， 收敛速度得到了明显的提高.为了进一步考察本文提出的改进AMDM的有效性，本文在非均匀Timoshenko梁的建模过程中考虑了横截面尺寸按线性函数、幂函数和指数函数变化的情形.通过数值算例检验方法的收敛性并且讨论影响收敛性的因素.此外，将计算结果与有限元仿真进行比较，验证本文方法的有效性和准确性.

1 模 型

图1为长度为l的弹性边界变截面梁，b_L和h_L分别为左侧边界的宽和高，b_R和h_R分别为右侧边界的宽和高，k_TL和k_TR分别为左侧和右侧边界的位移弹簧刚度，k_RL和k_RR分别为左侧和右侧边界的扭转弹簧刚度.定义形状函数b（x）和h（x）分别表示宽度和高度在轴向的变化，其截面面积和截面惯性矩可以表示为：

A(x)=A_{\mathrm{L}}{F}_{\mathrm{A}}(x)=b(x)h(x)（1）

I(x)=I_{\mathrm{L}}{F}_{\mathrm{I}}(x)=\frac{b(x)h^{\mathrm{3}}(x)}{\mathrm{12}}（2）

式中：A_L和I_L分别为左侧边界的截面面积和截面惯性矩； F_A（x）和F_I（x）分别表示位置x的面积和截面惯性矩相对于左侧边界（x=0）的截面面积A_L和截面惯性矩I_L的系数函数.

变截面Timoshenko梁的振动方程可表示为^［6］：

\frac{\partial }{\partial x}\left\{\kappa GA(x)\left[\frac{\partial y(x,t)}{\partial x}-\psi (x,t)\right]\right\}-\rho A(x)\frac{{\partial }^{\mathrm{2}}y(x,t)}{\partial t^{\mathrm{2}}}=\mathrm{0}（3）

\frac{\partial }{\partial x}\left[EI(x)\frac{\partial \psi (x,t)}{\partial x}\right]+\kappa GA(x)\left[\frac{\partial y(x,t)}{\partial x}-\psi (x,t)\right]-

\rho I(x)\frac{{\partial }^{\mathrm{2}}\psi (x,t)}{\partial t^{\mathrm{2}}}=\mathrm{0}（4）

式中：y(x,t)和\psi (x,t)分别为梁的横向位移和转角；A（x）和I（x）分别是梁的截面面积和惯性矩；\kappa为剪切系数；G为剪切模量；E为弹性模量；ρ为密度.

设\omega为振动的固有频率，Y(x)、\Psi (x)为模态位移和转角，利用y(x,t)=Y(x){\mathrm{e}}^{\mathrm{i}\omega t}和\psi (x,t)=\Psi (x){\mathrm{e}}^{\mathrm{i}\omega t}对式（3）和式（4）分离变量并进行无量纲化， 可得：

\frac{\mathrm{d}}{\mathrm{d}X}\left\{F_{\mathrm{A}}(X)\left[\frac{\mathrm{d}Y(X)}{\mathrm{d}X}-\Psi (X)\right]\right\}+\frac{{\Omega }^{\mathrm{2}}}{\xi }{F}_{\mathrm{A}}(X)Y(X)=\mathrm{0}（5）

\frac{\mathrm{d}}{\mathrm{d}X}\left[F_{\mathrm{I}}(X)\frac{\mathrm{d}\Psi (X)}{\mathrm{d}X}\right]+\xi F_{\mathrm{A}}(X)\left[\frac{\mathrm{d}Y(X)}{\mathrm{d}X}-\Psi (X)\right]+

\eta {\Omega }^{\mathrm{2}}{F}_{\mathrm{I}}(X)\Psi (X)=\mathrm{0}（6）

式中：X=\frac{x}{l}；Y(X)=\frac{Y(x)}{l}；\Psi (X)=\Psi (x)；F_{\mathrm{A}}(X)=F_{\mathrm{A}}(x)；F_{\mathrm{I}}(X)=F_{\mathrm{I}}(x)；{\Omega }^{\mathrm{2}}=\frac{\rho A_{\mathrm{L}}{\omega }^{\mathrm{2}}{l}^{\mathrm{4}}}{E{I}_{\mathrm{L}}}；\eta =\frac{I_{\mathrm{L}}}{A_{\mathrm{L}}{l}^{\mathrm{2}}}；\xi =\frac{\kappa G{A}_{\mathrm{L}}{l}^{\mathrm{2}}}{E{I}_{\mathrm{L}}}.\Omega为无量纲固有频率； \eta是与跨高比的平方成反比的无量纲参数； \xi是与剪切系数和跨高比相关的无量纲参数.

定义如下无量纲量：

\begin{array}{l}{K}_{\mathrm{T}\mathrm{L}}=\frac{k_{\mathrm{T}\mathrm{L}}{l}^{\mathrm{3}}}{E{I}_{\mathrm{L}}}\\ K_{\mathrm{R}\mathrm{L}}=\frac{k_{\mathrm{R}\mathrm{L}}l}{E{I}_{\mathrm{L}}}\\ K_{\mathrm{T}\mathrm{R}}=\frac{k_{\mathrm{T}\mathrm{R}}{l}^{\mathrm{3}}}{E{I}_{\mathrm{R}}}\\ K_{\mathrm{R}\mathrm{R}}=\frac{k_{\mathrm{R}\mathrm{R}}l}{E{I}_{\mathrm{R}}}\end{array}（7）

式中：K_TL、K_RL、K_TR和K_RR分别为左侧位移弹簧、左侧扭转弹簧、右侧位移弹簧和右侧扭转弹簧的无量纲刚度.

边界条件为：

\begin{array}{l}{\left.\left[\xi \frac{\mathrm{d}Y(X)}{\mathrm{d}X}-K_{\mathrm{T}\mathrm{L}}Y(X)-\xi \Psi (X)\right]\right|}_{X=\mathrm{0}}=\mathrm{0}\\ {\left.\left[\frac{\mathrm{d}\Psi (X)}{\mathrm{d}X}-K_{\mathrm{R}\mathrm{L}}\Psi (X)\right]\right|}_{X=\mathrm{0}}=\mathrm{0}\end{array}（8）

\begin{array}{l}{\left.\left[\begin{array}{l}\xi \frac{F_{\mathrm{A}}(X)}{F_{\mathrm{I}}(X)}\frac{\mathrm{d}Y(X)}{\mathrm{d}X}+K_{\mathrm{T}\mathrm{R}}Y(X)-\\ \xi \frac{F_{\mathrm{A}}(X)}{F_{\mathrm{I}}(X)}\Psi (X)\end{array}\right]\right|}_{X=\mathrm{1}}=\mathrm{0}\\ {\left.\left[\frac{\mathrm{d}\Psi (X)}{\mathrm{d}X}+K_{\mathrm{R}\mathrm{R}}\Psi (X)\right]\right|}_{X=\mathrm{1}}=\mathrm{0}\end{array}（9）

利用链式法则将式（5）和式（6）中的求导符号展开， 化简得：

\frac{{\mathrm{d}}^{\mathrm{2}}Y(X)}{\mathrm{d}{X}^{\mathrm{2}}}=-\frac{\mathrm{1}}{F_{\mathrm{A}}(X)}\frac{\mathrm{d}{F}_{\mathrm{A}}(X)}{\mathrm{d}X}\frac{\mathrm{d}Y(X)}{\mathrm{d}X}+

\frac{\mathrm{1}}{F_{\mathrm{A}}(X)}\frac{\mathrm{d}{F}_{\mathrm{A}}(X)}{\mathrm{d}X}\Psi (X)+\frac{\mathrm{d}\Psi (X)}{\mathrm{d}X}-\frac{{\Omega }^{\mathrm{2}}}{\xi }Y(X)（10）

\begin{array}{l}\frac{{\mathrm{d}}^{\mathrm{2}}\Psi (X)}{\mathrm{d}{X}^{\mathrm{2}}}=-\frac{\mathrm{1}}{F_{\mathrm{I}}(X)}\frac{\mathrm{d}{F}_{\mathrm{I}}(X)}{\mathrm{d}X}\frac{\mathrm{d}\Psi (X)}{\mathrm{d}X}-\\ \xi \frac{F_{\mathrm{A}}(X)}{F_{\mathrm{I}}(X)}\frac{\mathrm{d}Y(X)}{\mathrm{d}X}+\xi \frac{F_{\mathrm{A}}(X)}{F_{\mathrm{I}}(X)}\Psi (X)-\eta {\Omega }^{\mathrm{2}}\Psi (X)\end{array}（11）

将式（10）和式（11）中与形状函数有关的已知函数\frac{\mathrm{1}}{F_{\mathrm{A}}(X)}\frac{\mathrm{d}{F}_{\mathrm{A}}(X)}{\mathrm{d}X}，\frac{\mathrm{1}}{F_{\mathrm{I}}(X)}\frac{\mathrm{d}{F}_{\mathrm{I}}(X)}{\mathrm{d}X}和\frac{F_{\mathrm{A}}(X)}{F_{\mathrm{I}}(X)}定义为形状相关项， 其Taylor展开为：

\frac{\mathrm{1}}{F_{\mathrm{A}}(X)}\frac{\mathrm{d}{F}_{\mathrm{A}}(X)}{\mathrm{d}X}={\displaystyle \sum _{k=\mathrm{0}}^{\mathrm{\infty }}}{A}_k{X}^k

\frac{\mathrm{1}}{F_{\mathrm{I}}(X)}\frac{\mathrm{d}{F}_{\mathrm{I}}(X)}{\mathrm{d}X}={\displaystyle \sum _{k=\mathrm{0}}^{\mathrm{\infty }}}{B}_k{X}^k

\frac{F_{\mathrm{A}}(X)}{F_{\mathrm{I}}(X)}={\displaystyle \sum _{k=\mathrm{0}}^{\mathrm{\infty }}}{E}_k{X}^k（12）

式中：A_k 、B_k 、E_k 为X的k次项系数， 其表达式为：

{\left.A_k=\frac{\mathrm{1}}{k!}\frac{{\mathrm{d}}^k}{\mathrm{d}{X}^k}\left(\frac{\mathrm{1}}{F_{\mathrm{A}}(X)}\frac{\mathrm{d}{F}_{\mathrm{A}}(X)}{\mathrm{d}X}\right)\right|}_{X=\mathrm{0}}

B_k={\left.\frac{\mathrm{1}}{k!}\frac{{\mathrm{d}}^k}{\mathrm{d}{X}^k}\left(\frac{\mathrm{1}}{F_{\mathrm{I}}(X)}\frac{\mathrm{d}{F}_{\mathrm{I}}(X)}{\mathrm{d}X}\right)\right|}_{X=\mathrm{0}}

{\left.E_k=\frac{\mathrm{1}}{k!}\frac{{\mathrm{d}}^k}{\mathrm{d}{X}^k}\left(\frac{F_{\mathrm{A}}(X)}{F_{\mathrm{I}}(X)}\right)\right|}_{X=\mathrm{0}}（13）

利用Adomian分解法将横向位移和转角的模态函数Y(X)和\Psi (X)分解为级数形式^［25］：

Y(X)={\displaystyle \sum _{k=\mathrm{0}}^{\mathrm{\infty }}}{C}_k{X}^k

\Psi (X)={\displaystyle \sum _{k=\mathrm{0}}^{\mathrm{\infty }}}{D}_k{X}^k（14）

式中：C_k 和D_k 为X的k次项系数.

对式（14）求导可得：

\frac{\mathrm{d}Y(X)}{\mathrm{d}X}=\frac{\mathrm{d}}{\mathrm{d}X}{\displaystyle \sum _{k=\mathrm{0}}^{\mathrm{\infty }}}{C}_k{X}^k={\displaystyle \sum _{k=\mathrm{0}}^{\mathrm{\infty }}}(k+\mathrm{1})C_{k-\mathrm{1}}{X}^k（15）

\frac{\mathrm{d}\Psi (X)}{\mathrm{d}X}=\frac{\mathrm{d}}{\mathrm{d}X}{\displaystyle \sum _{k=\mathrm{0}}^{\mathrm{\infty }}}{D}_k{X}^k={\displaystyle \sum _{k=\mathrm{0}}^{\mathrm{\infty }}}(k+\mathrm{1})D_{k+\mathrm{1}}{X}^k（16）

将X=0代入式（14）~式（16）可得：

C_{\mathrm{0}}=Y(\mathrm{0}),C_{\mathrm{1}}=Y\text{'}(\mathrm{0})（17）

D_{\mathrm{0}}=\Psi (\mathrm{0}),D_{\mathrm{1}}=\Psi \text{'}(\mathrm{0})（18）

定义微分算符L={\mathrm{d}}^{\mathrm{2}}/\mathrm{d}{X}^{\mathrm{2}}， 则L的逆运算符L^-1为二重积分， L^{-\mathrm{1}}=\int \int (•)\mathrm{d}X\mathrm{d}X， 对式（10）和式（11）的等号两端应用L^-1， 并考虑式（12）~式（18）， 再利用级数的Cauchy乘积^［22］， 可得：

\begin{array}{l}Y(X)=C_{\mathrm{0}}+C_{\mathrm{1}}X+\\ L^{-\mathrm{1}}\left\{{\displaystyle \sum _{k=\mathrm{0}}^{\mathrm{\infty }}}{X}^k\left[\begin{array}{l}-{\displaystyle \sum _{j=\mathrm{0}}^k}(j+\mathrm{1})A_{k-j}{C}_{j+\mathrm{1}}+\\ {\displaystyle \sum _{j=\mathrm{0}}^k}{A}_j{D}_{k-j}+(k+\mathrm{1})D_{k+\mathrm{1}}-\\ \frac{{\Omega }^{\mathrm{2}}}{\xi }{C}_k\end{array}\right]\right\}\end{array}（19）

\begin{array}{l}\Psi (X)=D_{\mathrm{0}}+D_{\mathrm{1}}X+\\ L^{-\mathrm{1}}\left\{{\displaystyle \sum _{k=\mathrm{0}}^{\mathrm{\infty }}}{X}^k\left[\begin{array}{l}-{\displaystyle \sum _{j=\mathrm{0}}^k}(j+\mathrm{1})B_{k-j}{D}_{j+\mathrm{1}}-\\ \xi {\displaystyle \sum _{j=\mathrm{0}}^k}(j+\mathrm{1})E_{k-j}{C}_{j+\mathrm{1}}+\\ \xi {\displaystyle \sum _{j=\mathrm{0}}^k}{E}_j{D}_{k-j}-\eta {\Omega }^{\mathrm{2}}{D}_k\end{array}\right]\right\}\end{array}（20）

对式（19）~式（20）中的L^-1进行积分运算， 令m=k+2，比较式（14）与积分结果中X的同次幂项系数， 导出m≥2时，C_m 和D_m 的递推形式，其表达式为：

C_m=\frac{\mathrm{1}}{m(m-\mathrm{1})}\left[\begin{array}{l}-{\displaystyle \sum _{j=\mathrm{0}}^{m-\mathrm{2}}}(j+\mathrm{1})A_{m-\mathrm{2}-j}{C}_{j+\mathrm{1}}+\\ {\displaystyle \sum _{j=\mathrm{0}}^{m-\mathrm{2}}}{A}_j{D}_{m-\mathrm{2}-j}+\\ (m-\mathrm{1})D_{m-\mathrm{1}}-\frac{{\Omega }^{\mathrm{2}}}{\xi }{C}_{m-\mathrm{2}}\end{array}\right]（21）

D_m=\frac{\mathrm{1}}{m(m-\mathrm{1})}\left[\begin{array}{l}-{\displaystyle \sum _{j=\mathrm{0}}^{m-\mathrm{2}}}(j+\mathrm{1})B_{m-\mathrm{2}-j}{D}_{j+\mathrm{1}}-\\ \xi {\displaystyle \sum _{j=\mathrm{0}}^{m-\mathrm{2}}}(j+\mathrm{1})E_{m-\mathrm{2}-j}{C}_{j+\mathrm{1}}+\\ \xi {\displaystyle \sum _{j=\mathrm{0}}^{m-\mathrm{2}}}{E}_j{D}_{m-\mathrm{2}-j}-\eta {\Omega }^{\mathrm{2}}{D}_{m-\mathrm{2}}\end{array}\right]（22）

由于无法得到无穷级数的精确解，因此，需要对式（14）进行截断近似^［23］. 令Y^{[n]}(X)={\displaystyle \sum _{k=\mathrm{0}}^{n-\mathrm{1}}}{C}_k{X}^k， {\Psi }^{[n]}(X)={\displaystyle \sum _{k=\mathrm{0}}^{n-\mathrm{1}}}{D}_k{X}^k表示模态函数的n阶近似解.相应地，形状相关项式（12）的Taylor展开式也取前n项.式（19）和式（20）中的初始系数可以由边界条件确定.通过观察边界条件表达式（8）和式（9）的取值，并且根据不同的边界条件设置任意常数a₁和a₂，可以得到弹性边界变截面Timoshenko梁的边界条件， 见表1和表2.

利用表1确定初始系数C₀、C₁、D₀和D₁， 代入递推式（21）和式（22）可以依次求出剩余的C_m 和D_m， 从而确定模态函数. 将模态函数代入表2中的边界条件， 提取a₁和a₂的系数， 可得：

\left[\begin{array}{cc}{F}_{\mathrm{11}}\left({\Omega }^{\mathrm{2}}\right)& F_{\mathrm{12}}\left({\Omega }^{\mathrm{2}}\right)\\ F_{\mathrm{21}}\left({\Omega }^{\mathrm{2}}\right)& F_{\mathrm{22}}\left({\Omega }^{\mathrm{2}}\right)\end{array}\right]\left[\begin{array}{l}{a}_{\mathrm{1}}\\ a_{\mathrm{2}}\end{array}\right]=\mathrm{0}（23）

式中：F₁₁、F₁₂、F₂₁、F₂₂都是Ω²的函数，根据Cramer法则^［22］，式（23）的系数行列式必为零，由此可导出特征方程.求解特征方程能够得到第i阶无量纲固有频率Ω_i .

为了解决传统AMDM求解变截面Timoshenko梁的振动问题时，在整个梁的跨度内收敛性差的问题，本文对梁进行分段，并且定义基于IAE的误差指标：

F_{\mathrm{E}\mathrm{r}\mathrm{r}}^{[n]}={\int }_{\mathrm{0}}^{L_{\mathrm{s}\mathrm{e}\mathrm{g}}}\left|F_{\mathrm{T}\mathrm{a}\mathrm{y}\mathrm{l}\mathrm{o}\mathrm{r}}^{[n]}(X)-F_{\mathrm{O}\mathrm{r}\mathrm{g}}(X)\right|\mathrm{d}X（24）

式中：n为Taylor展开式的项数；F_{\mathrm{O}\mathrm{r}\mathrm{g}}是原函数；F_{\mathrm{T}\mathrm{a}\mathrm{y}\mathrm{l}\mathrm{o}\mathrm{r}}^{[n]}是保留n项的Taylor展开式；L_{\mathrm{s}\mathrm{e}\mathrm{g}}是每段无量纲长度.通过逐渐缩小分段长度使IAE满足特定阈值，以便能够根据变截面梁的锥度和幂次自动确定分段的长度和数量.

在计算每个分段时，需要进行坐标平移， 使得左端点坐标始终为0，右端点坐标为第j个分段的长度L_j .左侧端点的高度和宽度h₀和b₀要相应地修改为h_j-1和b_j-1.

2 形状函数与形状相关项

本节选取常见的几种变截面梁的形状函数进行分析.

2.1 线性形状函数

线性变截面梁的形状函数为^［23］：

b(X)=b_{\mathrm{L}}(\mathrm{1}-{\beta }_{\mathrm{b}}X)（25）

h(X)=h_{\mathrm{L}}(\mathrm{1}-{\beta }_{\mathrm{h}}X)（26）

式中：β_b=1-b_R/b_L；β_h=1-h_R/h_L；b_L和b_R分别为梁的左、右侧边界的宽度；h_L和h_R分别为梁的左、右侧边界的高度. β_b和β_h分别为梁沿宽、高方向的锥度， β_b和β_h越大， 锥度越大. 截面面积和惯性矩系数函数为：

F_{\mathrm{A}}(X)=\left(\mathrm{1}-{\beta }_{\mathrm{b}}X\right)\left(\mathrm{1}-{\beta }_{\mathrm{h}}X\right)（27）

F_{\mathrm{I}}(X)=\left(\mathrm{1}-{\beta }_{\mathrm{b}}X\right){\left(\mathrm{1}-{\beta }_{\mathrm{h}}X\right)}^{\mathrm{3}}（28）

对应的形状相关项Taylor展开式系数为：

A_k=-{\beta }_{\mathrm{b}}^{k+\mathrm{1}}-{\beta }_{\mathrm{h}}^{k+\mathrm{1}}B_k=-{\beta }_{\mathrm{b}}^{k+\mathrm{1}}-\mathrm{3}{\beta }_{\mathrm{h}}^{k+\mathrm{1}}（29）

E_k=(k+\mathrm{1}){\beta }_{\mathrm{h}}^k

2.2 幂函数变截面梁的形状函数

幂函数变截面梁的形状函数为^［2-3］：

h(X)=h_{\mathrm{L}}(\mathrm{1}-{\beta }_{\mathrm{h}}{X}^p)（30）

b(X)=b_{\mathrm{L}}(\mathrm{1}-{\beta }_{\mathrm{b}}{X}^p)（31）

式中： p为变截面梁函数的幂次. 截面面积和惯性矩系数函数为：

F_{\mathrm{A}}(X)=(\mathrm{1}-{\beta }_{\mathrm{b}}{X}^p)(\mathrm{1}-{\beta }_{\mathrm{h}}{X}^p)（32）

F_{\mathrm{I}}(X)=(\mathrm{1}-{\beta }_{\mathrm{b}}{X}^p)(\mathrm{1}-{\beta }_{\mathrm{h}}{X}^p{)}^{\mathrm{3}}（33）

对应的形状相关项Taylor展开式系数为：

A_k=\left[\mathrm{1}-{\left(-\mathrm{1}\right)}^k\right]\left(\begin{array}{l}\mathrm{2}{\displaystyle \sum _{j=\mathrm{1}}^{k/\mathrm{2}-\mathrm{1}}}{\beta }_{\mathrm{b}}^j{\beta }_{\mathrm{h}}^{k/\mathrm{2}-j}-\\ \left({\beta }_{\mathrm{b}}+{\beta }_{\mathrm{h}}\right){\displaystyle \sum _{j=\mathrm{0}}^{k/\mathrm{2}-\mathrm{1}}}{\beta }_{\mathrm{b}}^j{\beta }_{\mathrm{h}}^{k/\mathrm{2}-\mathrm{1}-j}\end{array}\right)

B_k=\left[\mathrm{1}-{\left(-\mathrm{1}\right)}^k\right]\left(\begin{array}{l}\mathrm{4}{\displaystyle \sum _{j=\mathrm{1}}^{k/\mathrm{2}-\mathrm{1}}}{\beta }_{\mathrm{b}}^j{\beta }_{\mathrm{h}}^{k/\mathrm{2}-j}-\\ \left({\beta }_{\mathrm{b}}+\mathrm{3}{\beta }_{\mathrm{h}}\right){\displaystyle \sum _{j=\mathrm{0}}^{k/\mathrm{2}-\mathrm{1}}}{\beta }_{\mathrm{b}}^j{\beta }_{\mathrm{h}}^{k/\mathrm{2}-\mathrm{1}-j}\end{array}\right)

E_k=\frac{\mathrm{1}}{\mathrm{2}}\left[\mathrm{1}+{\left(-\mathrm{1}\right)}^k\right]\left(\frac{k}{\mathrm{2}}+\mathrm{1}\right){\beta }_{\mathrm{h}}^{k/\mathrm{2}}（34）

2.3 指数函数变截面梁的形状函数

指数函数变截面梁的形状函数为^［2］：

h(X)=h_{\mathrm{L}}{\mathrm{e}}^{{\delta }_{\mathrm{h}}X}（35）

b(X)=b_{\mathrm{L}}{\mathrm{e}}^{{\delta }_{\mathrm{b}}X}（36）

式中：δ_b和δ_h分别为梁的宽、高方向的锥度. 截面面积和惯性矩系数函数为：

F_{\mathrm{A}}(X)={\mathrm{e}}^{\left({\delta }_{\mathrm{b}}+{\delta }_{\mathrm{h}}\right)X}（37）

F_{\mathrm{I}}(X)={\mathrm{e}}^{\left({\delta }_{\mathrm{b}}+\mathrm{3}{\delta }_{\mathrm{h}}\right)X}

对应的形状相关项Taylor展开式系数为：

A_{\mathrm{0}}={\delta }_{\mathrm{b}}+{\delta }_{\mathrm{h}}（38）

B_{\mathrm{0}}={\delta }_{\mathrm{b}}+\mathrm{3}{\delta }_{\mathrm{h}}E_k={(-\mathrm{2})}^k\frac{{\delta }_{\mathrm{h}}^k}{k!}

3 数值算例

对于各向同性材料， 其剪切模量G与弹性模量E满足：

G=\frac{E}{\mathrm{2}(\mathrm{1}+v)}（39）

式中：v为泊松比.

本文所有算例均采用v=0.25， \kappa =\mathrm{5}/\mathrm{6}. 计算使用的CPU为I9-9900KF， 内存为16 G.

3.1 自适应分段的效果

在进行计算之前， 需要根据形状相关项的收敛情况进行自适应分段， 为了显示自适应分段的有效性， 图2给出幂函数变截面梁的形状相关项Taylor展开式自适应分段与不分段的对比结果， 自适应分段点已在图中标出. 由图2可见， 不分段时形状相关项的Taylor展开式与原函数偏差较大， 经过自适应分段的Taylor展开式与原函数完全吻合. 线性和指数函数变截面梁的计算结果与此相似.

图3给出变截面梁前5阶固有频率收敛精度为0.000 1时的计算耗时.图3（a）给出3种变截面梁的计算耗时随锥度的变化情况.在不分段的情况下， 计算耗时都随锥度增加而增加，其中，幂函数变截面梁计算耗时上升最快，指数函数变截面梁计算耗时上升趋势最为平缓.采用自适应分段的计算耗时明显减少，在锥度为0.8的情况下，线性、幂函数和指数函数变截面梁的计算耗时分别缩短了69%、91%和46%.图3（b）给出锥度0.5的幂函数变截面梁的计算耗时随幂次的变化情况.结果表明， 幂次越高，计算耗时越长，自适应分段能够明显地缩短计算时间，特别是在幂次为5时，与不分段相比计算耗时减少了70%.图3（c）表示锥度0.7的幂函数变截面梁自适应分段与平均分段计算耗时的对比结果， 由此可见，虽然对于低幂次形状函数， 均分为5段的收敛性能与自适应分段结果相近， 但平均分段无法预估分段长度的合理性， 而自适应分段能够根据幂次、锥度合理地确定分段长度和数量， 从而提高计算效率， 特别是当幂次高达5时， 自适应分段的计算耗时比均分5段缩短了38%， 比均分10段缩短了65%.

3.2 方法的收敛性

表3表示幂函数（p=2）变截面梁的Adomian分解最高幂次n取不同值时，前5阶无量纲固有频率的收敛情况，各阶固有频率的收敛数值已被加粗.计算过程基本参数：K_TL=K_RL=K_TR=K_RR=100， β_b=β_h=0.5，l/h=10.从表3可见，3~5阶固有频率随着n的增大先后收敛，低阶固有频率的收敛速度相对更快.当n为38时，前5阶固有频率均收敛至预设精度（0.000 01）.线性和指数函数变截面梁的收敛情况与此相似.

3.3 方法的准确性

表4给出幂函数（p=2）变截面梁前5阶无量纲固有频率与有限元仿真结果的比较，输入基本参数为： K_TL=K_RL=K_TR=K_RR， β_b=β_h=0.5，l/h=10.能够看到， 表4中的相对误差均不超过0.16%.线性和指数函数变截面梁的计算结果与此相似.此外，在两端自由的情况下， 理论上会出现前两阶固有频率均为零的刚体模态. 有限元仿真得到的固有频率只能近似为零， 无法达到完全为零的精确结果， 而本文的方法可精确预测梁的刚体模态频率.

3.4 弹簧类型对固有频率的影响

图4表示线性、幂函数和指数函数变截面梁弹簧类型对固有频率的影响（l/h=10， K_TL=K_RL=100， 锥度0.5），纵坐标定义为固有频率与右侧自由梁固有频率的比值.由图4可见，随着刚度K_TR和K_RR的增加，各阶固有频率相对于右侧自由梁固有频率的比值曲线单调上升；只施加位移弹簧K_TR能够获得比只施加扭转弹簧K_RR更高的固有频率比极值和更陡峭的上升曲线，即位移弹簧对固有频率的影响强于扭转弹簧；同时施加位移弹簧和扭转弹簧，固有频率随弹簧刚度的上升曲线比只施加一种弹簧更加陡峭， 达到的极值也更高.弹簧刚度对低阶固有频率的影响强于对高阶固有频率的影响.幂函数变截面梁的固有频率受边界条件影响最小，指数函数变截面梁的固有频率受边界条件影响最大.另外，各阶固有频率均随弹簧刚度的增加而增加，变化趋势为先缓后急， 随后趋于平缓， 最终趋近于一端固定梁的固有频率.

3.5 弹簧刚度对固有频率-锥度曲线的影响

图5表示锥度对弱弹性边界条件变截面梁固有频率的影响，其中，纵坐标定义为固有频率与等截面梁固有频率的比值.由图5可见， 随锥度的增加，第1阶和第2阶固有频率呈现先减后增的趋势， 曲线拐点在β_h=0.4附近；幂函数变截面梁的第3阶固有频率也呈现先减后增的趋势， 但幂次越高拐点越不明显， 线性和指数函数变截面梁的第3阶固有频率有单调减小的趋势.

图6表示锥度对较强弹性边界条件变截面梁固有频率的影响. 由图6可见，当弹簧刚度增加时，各阶固有频率随锥度的增加有总体下降的趋势，在锥度较大（β_h>0.7）时，除了指数函数变截面梁外，各阶固有频率均有向上的拐点.

3.6 跨高比对固有频率的影响

图7表示跨高比对不同阶次固有频率的影响. 取K_TL=K_RL=K_TR=K_RR=100， l/h=10， 对于矩形梁， 无量纲参数η=（h/l）²/12，与跨高比的平方成反比.图7纵坐标为无量纲固有频率与Euler-Bernoulli梁无量纲固有频率Ω_i（EB）的比值.由图7可见，随着η的增加（跨高比减小），固有频率减小，当η较小时，Timo-shenko梁与Euler-Bernoulli梁的固有频率趋于一致， 即剪切变形和转动惯量可以忽略.幂函数变截面梁受跨高比的影响最大，指数函数变截面梁受跨高比的影响最小，这是因为当锥度相等时，幂函数变截面梁最粗，指数函数变截面梁最细.粗梁受剪切变形的影响较大，而细梁受剪切变形的影响较小.

4 结 论

本文采用自适应分段技术对AMDM进行改进， 得到了一种求解具有不同弹性边界的变截面Timoshenko梁振动问题近似解的有效方法. 将改进的AMDM应用于幂函数和指数函数形式的复杂变截面Timoshenko梁， 仿真结果表明， 极大地改善和优化了近似求解的收敛性能， 计算耗时在不同情形下缩短了46%~91%. 采用改进的AMDM得到的固有频率近似解与有限元仿真结果相比， 误差可控制在不超过0.16%.

1）当位移和扭转弹簧刚度较小时， 低阶固有频率随锥度的增大而增大， 高阶固有频率随锥度的增大而降低， 当位移和扭转弹簧刚度较大时， 各阶固有频率都会随着锥度的增大而降低.

2）跨高比的增加会导致剪切变形的影响减小， 因此， Timoshenko梁的固有频率随之逐步增加， 并最终趋向于Euler-Bernoulli梁的固有频率.

3）跨高比对幂函数变截面梁影响最大， 对指数函数变截面梁的影响最小.

参考文献

[1] 马静敏， 任勇生， 姚文莉 ．复合材料变截面旋转薄壁悬臂梁自由振动分析［J］．工程力学， 2013， 30（1）：37-44．

[2] MA J M， REN Y S， YAO W L. Free vibration analysis of rotating composite thin-walled cantilever beams with variable closed-section［J］. Engineering Mechanics， 2013， 30（1）：37-44．（in Chinese）

[3] MAHMOUD M A. Natural frequency of axially functionally graded， tapered cantilever beams with tip masses［J］. Engineering Structures， 2019， 187：34-42．

[4] 张永康， 鲍四元 ．幂函数形式连续变截面梁振动的弯曲固有频率分析［J］．应用声学， 2024， 43（2）：330-338．

[5] ZHANG Y K， BAO S Y ．Analysis of natural frequency of beams with power-function-type continuous variable section［J］. Journal of Applied Acoustics， 2024， 43（2）： 330-338．（in Chinese）

[6] 周坤涛， 杨涛， 葛根， 等 ．基于Bessel和Meijer-G函数的楔形和锥形悬臂梁振动分析［J］．振动与冲击， 2022， 41（4）：253-261．

[7] ZHOU K T， YANG T， GE G， et al ．Vibration analysis of wedge and cone cantilever beams based on Bessel and Meijer-G functions［J］. Journal of Vibration and Shock， 2022， 41（4）： 253-261．（in Chinese）

[8] 牛国华， 王刚锋， 王剑 ．组合L型变截面梁的自由振动特性分析［J］．振动与冲击， 2022， 41（5）： 228-234．

[9] NIU G H， WANG G F， WANG J ．Analysis of free vibration characteristics of composite L-shaped variable cross-section beam［J］．Journal of Vibration and Shock， 2022， 41（5）：228-234．（in Chinese）

[10] TIMOSHENKO S， YOUNG D H， WEAVER W. Vibration problems in engineering［EB/OL］.（2016-07-04）［2025-03-17］.

[11] 潘旦光， 吴顺川， 张维 ．变截面Timoshenko悬臂梁自由振动分析［J］．土木建筑与环境工程， 2009， 31（3）： 25-28．

[12] PAN D G， WU S C， ZHANG W. Free vibration of non-uniform Timoshenko cantilever beams［J］. Journal of Civil， Architectural & Environmental Engineering， 2009， 31（3）：25-28．（in Chinese）

[13] 葛仁余， 张金轮， 姜忠宇， 等. 轴向功能梯度变截面Timoshenko梁自由振动的研究［J］. 振动与冲击， 2017， 36（22）： 158-165．

[14] GE R Y， ZHANG J L， JIANG Z Y， et al ．Free vibration analysis of axially functionally Timoshenko beams with a non-uniform cross-section［J］．Journal of Vibration and Shock， 2017， 36（22）：158-165．（in Chinese）

[15] 黄意新， 穆洲， 郭明全， 等 ．复杂边界条件轴向功能梯度梁动力学分析［J］. 哈尔滨工业大学学报， 2018， 50（10）： 143-150．

[16] HUANG Y X， MU Z， GUO M Q， et al ．Dynamic analysis of axially functionally graded beams with complex boundary conditions［J］. Journal of Harbin Institute of Technology， 2018， 50（10）： 143-150．（in Chinese）

[17] MOLINA-VILLEGAS J C， BALLESTEROS ORTEGA J E， SOTO S B ．Closed-form solutions for axially non-uniform Timoshenko beams and frames under static loading［J］. Composite Structures， 2024， 337：118078．

[18] KARAMI G， MALEKZADEH P， SHAHPARI S A. A DQEM for vibration of shear deformable nonuniform beams with general boundary conditions［J］. Engineering Structures， 2003， 25（9）：1169-1178．

[19] AUCIELLO N M. Free vibration of a restrained shear-deformable tapered beam with a tip mass at its free end［J］. Journal of Sound and Vibration， 2000， 237（3）：542-549.

[20] AUCIELLO N M， ERCOLANO A. A general solution for dynamic response of axially loaded non-uniform Timoshenko beams［J］．International Journal of Solids and Structures， 2004， 41（18/19）：4861-4874．

[21] GOEL R P ．Transverse vibrations of tapered beams［J］. Journal of Sound and Vibration， 1976， 47（1）：1-7．

[22] GROSSI R O， ARANDA A， BHAT R B ．Vibration of tapered beams with one end spring hinged and the other end with tip mass［J］．Journal of Sound and Vibration， 1993， 160（1）：175-178．

[23] GROSSI R O， DEL V ARENAS B ．A variational approach to the vibration of tapered beams with elastically restrained ends［J］. Journal of Sound and Vibration， 1996， 195（3）： 507-511．

[24] HO S H， CHEN C K. Analysis of general elastically end restrained non-uniform beams using differential transform［J］. Applied Mathematical Modelling， 1998， 22（4/5）： 219-234．

[25] DE ROSA M A， AUCIELLO N M. Free vibrations of tapered beams with flexible ends［J］. Computers & Structures， 1996， 60（2）： 197-202．

[26] AUCIELLO N M ．A comment on “a note on vibrating tapered beams” ［J］. Journal of Sound and Vibration， 1995， 187（4）：724-726．

[27] ADOMIAN G ．A review of the decomposition method in applied mathematics［J］. Journal of Mathematical Analysis and Applications， 1988， 135（2）： 501-544．

[28] ADOMIAN G ．A review of the decomposition method and some recent results for nonlinear equations［J］. Computers & Mathematics with Applications， 1991， 21（5）： 101-127．

[29] LAI H， CHEN C K， HSU J ．Free vibration of non-uniform Euler-bernoulli beams by the adomian modified decomposition method［J］. Computer Modeling in Engineering & Sciences， 2008， 34（1）： 87-116．

[30] ZHANG C J， REN Y S， JI S J， et al ．Analysis of the vibration and chatter stability of a tapered composite boring bar［J］．Shock and Vibration， 2020， 2020（1）：4190806．

[31] ZHOU H Y， HUANG J C， FU Y F ．Transverse free vibration of non-uniform beams with exponentially varying cross-sections［J］．International Journal of Applied Mechanics， 2025， 17（2）：2550010．

[32] HSU J C， LAI H Y， CHEN C K ．An innovative eigenvalue problem solver for free vibration of uniform Timoshenko beams by using the Adomian modified decomposition method［J］．Journal of Sound and Vibration， 2009， 325（1/2）：451-470．

[33] 郭银凤 ．改进的Adomian分解法在非线性方程中的应用［D］．大连：大连理工大学， 2019．

[34] GUO Y F ．Applications of modified adomian decomposition method to nonlinear equations［D］. Dalian：Dalian University of Technology， 2019．（in Chinese）

[35] GHOSH S， ROY A， ROY D. An adaptation of adomian decomposition for numeric-analytic integration of strongly nonlinear and chaotic oscillators［J］. Computer Methods in Applied Mechanics and Engineering， 2007， 196：1133-1153．

[36] JIAO Y C， YAMAMOTO Y， DANG C， et al ．An aftertreatment technique for improving the accuracy of Adomian’s decomposition method［J］．Computers & Mathematics with Applications， 2002， 43（6/7）： 783-798．

[37] 王永亮 ．变截面变曲率梁振型的有限元超收敛拼片恢复解和网格自适应分析［J］．工程力学， 2020， 37（12）：1-8．

[38] WANG Y L ．Superconvergent patch recovery solutions and adaptive mesh refinement analysis of finite element method for the vibration modes of non-uniform and variable curvature beams［J］．Engineering Mechanics， 2020， 37（12）： 1-8．（in Chinese）

[39] 徐亚飞， 肖映雄， 吴宇航 ．基于Python-Abaqus的自适应网格重划分算法实现及其应用［J］．计算力学学报， 2023， 40（5）：807-814．

[40] XU Y F， XIAO Y X， WU Y H. Implementation and its application of adaptive remeshing algorithm based on Python-Abaqus［J］. Chinese Journal of Computational Mechanics， 2023， 40（5）： 807-814．（in Chinese）

[41] PIETRALA M， LEŚNIEWSKI P， BARTOSZEWICZ A ．IAE minimization in sliding mode control with input and velocity constraints［J］．IEEE Access， 2022， 10：28631-28641．

基金资助

国家自然科学基金资助项目(11672166)

National Natural Science