功能梯度梁临界荷载与固有频率的微分求积法计算

熊海超，葛仁余，马国强，刘小双，余本源

(安徽工程大学 建筑工程学院，安徽 芜湖 241000)

功能梯度材料是一种新型材料，它是将多种性能各异的材料按照设计意愿形成性能连续变化的合成材料，以满足各种特殊工程结构的需要。但是，由于功能梯度材料的物理性质连续变化，给其力学分析带来一定的难题，除非一些特殊情况，一般用解析方法分析功能梯度材料的振动问题十分困难，如文献[1]用能量法获得了锥形变截面桩屈曲临界荷载的解析解。目前，数值方法是用来研究轴向功能梯度变截面梁的自由振动和屈曲临界荷载的主要手段，文献[2-3]运用精确的动态刚度方法，获得了Euler-Bernoulli梁和Timoshenko梁的固有频率和屈曲临界荷载；文献[4]运用微分求积单元法研究了轴向功能梯度Euler-Bernoulli梁自由振动和屈曲临界荷载的问题；基于Euler-Bernoulli梁理论和Timoshenko梁理论，文献[5-6]运用插值矩阵法研究了轴向功能梯度变截面梁自由振动和屈曲临界荷载问题。将边界条件转换为简便的形式，避免了刚度系数矩阵的无穷大值，文献[7]研究了轴向功能梯度变截面弹性地基梁的屈曲临界荷载问题。文献[8]采用局部微分求积法研究了带有弹性约束的轴向功能梯度变截面桩屈曲临界荷载问题。利用均匀梁单元的形状函数，提出了一种新的梁单元。文献[9]采用有限元计算方法，研究了轴向功能梯度梁的自由振动和稳定性问题。

研究采用切比雪夫多项式的根作为非均匀网格布点方式，由微分求积法(DQM)研究受轴向荷载作用的功能梯度材料Euler-Bernoulli变截面梁的屈曲临界荷载、固有频率，以及轴向荷载对固有频率的影响。基于Euler-Bernoulli梁理论，获得求解功能梯度材料Euler-Bernoulli梁的屈曲临界荷载和固有频率的变系数常微分方程，再根据微分求积法理论，将变系数常微分方程的特征值问题转化为线性代数方程组的特征值问题，可一次性地获得梁的屈曲临界荷载、固有频率，以及轴向荷载对梁的第1阶固有频率的影响。

1 计算模型

1.1 Euler-Bernoulli梁稳定性的基本理论

图1 轴向功能梯度变截面Euler-Bernoulli梁

轴向功能梯度变截面Euler-Bernoulli梁如图1所示。功能梯度材料Euler-Bernoulli变截面梁的长度为

l

，材料性能和截面面积沿轴向

x

任意连续变化。设梁自由振动时的挠度为

w(x,t)

，材料的弹性模量为

E(x)=E

f

(x)

，材料的密度为

ρ(x)=ρ

f

(x)

，截面面积为

A(x)=A

f

(x)

，截面转动惯量为

I(x)=I

f

(x)

，它们均为

x

的连续函数。其中，

E

、

ρ

、

A

、

I

对应于变截面梁在左端边界

x=

0位置处材料的弹性模量、密度、截面面积和截面惯性矩；

f

(x)

为

x

的函数

(i=

1

,

2

,

3

,

4

)

；设轴向荷载

P>

0为压力，

P<

0为拉力，则在轴向荷载

P

作用下的功能梯度材料Euler-Bernoulli变截面梁的自由振动控制方程为9：

(

1

)

研究主要考虑梁的谐波振动，则有：

w(x,t)=W(x)

sin

ωt,

(

2

)

式中，

ω

为梁振动角频率

;t

为时间变量

;w(x,t)

为梁挠度

;W(x)

为梁振型函数。联立式

(

1

)

和式

(

2

)

，可得轴向功能梯度Euler-Bernoulli变截面梁在轴向荷载

P

作用下的振动控制方程：

(

3

)

(

4

)

将式

(

4

)

代入式

(

3

)

，可得轴向荷载

P

作用下轴向功能梯度Euler-Bernoulli变截面梁的自由振动控制方程：

(

5

)

当轴向荷载

P

达到屈曲临界荷载时，功能梯度材料Euler-Bernoulli变截面梁的固有频率

ω=

0，将式

(

5

)

转化为式

(

6

)

关于屈曲临界荷载

λ

的变系数常微分方程的特征值问题：

(

6

)

1.2 Euler-Bernoulli梁的边界条件

Euler-Bernoulli梁的转角

θ

，弯矩

M

和剪力

T

分别为：

(

7

)

简支-简支梁

(

H-H

)

的边界条件是两端挠度和弯矩为零，即

W(

0

)=

0

,W

(

0

)=

0

;W(

1

)=

0

,W

(

1

)=

0

,

(8a)

固端-固端梁

(

C-C

)

的边界条件是两端挠度和转角为零，即

W(

0

)=

0

,W

(

0

)=

0

;W(

1

)=

0

,W

(

1

)=

0

,

(8b)

固端-简支梁

(

C-H

)

的边界条件是固端挠度和转角为零，简支端挠度和弯矩为零，即

W(

0

)=

0

,W

(

0

)=

0

;W(

1

)=

0

,W

(

1

)=

0

,

(8c)

因此，功能梯度材料Euler-Bernoulli变截面梁固有频率的计算实际成为求解边界条件式

(

8

)

下变系数常微分方程式

(

5

)

的特征值问题；屈曲临界荷载的计算实际成为求解边界条件式

(

8

)

下变系数常微分方程式

(

6

)

的特征值问题。研究采用微分求积法对上述问题展开数值方法分析。

2 微分求积法数值计算方法

功能梯度材料Euler-Bernoulli变截面梁微分求积法计算模型如图2所示。研究采用切比雪夫多项式的根作为节点坐标的方式如式(9a)所示，均匀等步长分布节点坐标的方式如式(9b)所示：

图2 轴向功能梯度材料Euler-Bernoulli梁的网格划分模型

(

9a

)

(

9b

)

考虑一维函数

W(ξ)

在区间

[

0

,

1

]

上可微，将区间

[

0

,

1

]

划分为

n

段，0

=ξ

，

ξ

，

ξ

，…

,ξ

-1，

ξ

=

1，在每个节点上函数及其导数用区间

[

0

,

1

]

上

n

+1个节点函数值的加权线性和近似表示。这里将梁的振型函数

W(ξ)

用拉格朗日

(

Lagrange

)

插值函数来描述，即

(

10

)

式中，

l

(ξ)

为拉格朗日多项式，其表达式为：

(

11

)

由式

(

10

)

对函数

W(ξ)

求一阶导数，得到：

(

12

)

将式

(

12

)

在

n

段区间

[

0

,

1

]

上进行离散，从而进一步得到：

(

13

)

依次类推可得：

(

14

)

由式

(

14

)

可知：

(

15

)

(

16

)

(

17

)

由于

(

18

)

因此，各阶导数的加权系数矩阵之间的关系为：

(

19

)

将式

(

6

)

计算屈曲临界荷载的微分方程组中变系数写成对角阵向量形式，

I

为单位阵，即

(

20

)

将式

(

6

)

微分方程组写成代数方程组向量形式为：

C

W

-

λC

W=

0

,

(

21

)

不失一般性，以两端固支梁

(

C-C

)

情况为例进行讨论，则相应的边界条件向量形式可表示为：

(

22

)

这里，

I

是0和

n

，

[

…

]

为矩阵

[

…

]

的第

I

行元素。研究求解如式

(

6

)

的四阶微分方程时，运用节点替代法

(δ

法

)

来处理边界条件，即将边界条件利用与梁两端相隔

δ

的两个点来替代

(

10

<δ<

10

)

，取节点序列为

{ξ}={

0

,δ,

…

,ξ

,

…

,

1-

δ,

1

}

(j=

3

,

4

,

…

,n

-2

),

(

23

)

(

24

)

3 数值算例和问题讨论

设轴向功能梯度Euler-Bernoulli变截面梁的长度为

l

，截面形状为矩形，材料的弹性模量和密度分别为：

梁的截面面积和惯性矩分别为：

式中，

c

和

c

分别是梁截面的宽、高锥度系数

;x

是以梁的左端为起点沿轴线方向的横坐标；当

c

=c

=

0时为等截面梁。运用微分求积法计算轴向功能梯度Euler-Bernoulli变截面梁屈曲临界荷载的变系数常微分方程式，即式

(

6

)

，离散单元数取

n=

20，非均匀网格

(

切比雪夫多项式根作为节点

)

和均匀网格两种布点方式的微分求积法计算值对比如表1所示。由表1计算结果可知，切比雪夫多项式根作为布点方式的计算值与文献

[

9

]

结果吻合良好，而均匀等步长布点方式的计算值与文献

[

9

]

结果误差较大甚至失真。大量的数值实验研究表明，微分求积法求解变系数微分方程时，采用式(9

b

)均匀节点划分网格时，微分求积法计算数值不稳定甚至失真，而一些非均匀节点划分网格往往能得到很好的结果，如切比雪夫多项式的根作为离散节点分布形式(式(9

a

))。

表1 在不同布点方式下，轴向功能梯度变截面梁的无量纲屈曲临界荷载计算值比对

在C-C、C-H与H-H不同边界条件和梁截面锥度系数下，以切比雪夫多项式根作为变步长非均匀网格布点方式，取

n=

20，微分求积法求解变系数常微分方程式

(

5

)

获得的荷载

λ

与第1阶固有频率

μ

的关系曲线图如图3所示。由图3可知：①随着轴向荷载由负值向正值增大，即轴向荷载由拉力向压力转化，梁的第1阶固有频率逐渐减小；当轴向荷载较小时，梁振动频率近似呈线性降低，当轴向荷载即将接近屈曲临界荷载值时，梁的第1阶横向振动频率开始非线性急剧降低；②当边界条件为H-H时，梁的边界约束最弱，轴向荷载对第1阶固有频率的影响明显，屈曲临界荷载最小；当边界条件为C-C时，梁的边界约束最强，屈曲临界荷载最大；③随着压力的增加，梁第1阶固有频率逐渐趋于零，当第1阶固有频率为零时，对应的轴向荷载即为梁的屈曲临界荷载，在边界条件C-C、C-H和H-H下，截面宽锥度系数

c

=

0.2和截面高锥度系数

c

=

0.2时，梁的无量纲屈曲临界荷载

λ

分别为37.602 327 3、19.335 410 3、9.597 120 84；截面宽锥度系数

c

=

0.8和截面高锥度系数

c

=

0.8时，无量纲屈曲临界荷载

λ

分别为21.781 289 8、11.218 322 3、5.622 775 59；截面宽锥度系数

c

=0.4和截面高锥度系数

c

=

0.4时，无量纲屈曲临界荷载

λ

分别为0.664 879 02、1.371 282 41、0.707 535 808，微分求积法求解变系数常微分方程式

(

5

)

的计算值与表1的计算值自行吻合。

表2 轴向功能梯度材料Euler-Bernoulli变截面梁的无量纲屈曲临界荷载计算值

图3 轴向荷载对Euler-Bernoulli梁第1阶固有频率的影响示意图

4 结论

研究基于Euler-Bernoulli梁理论，建立了求解功能梯度材料变截面梁的屈曲临界荷载和固有频率的变系数常微分方程，再运用微分求积法理论，通过缜密的数学推导，将变系数常微分方程的特征值问题转化为线性代数方程组的特征值问题，一次性地获得梁的屈曲临界荷载以及轴向力与梁固有频率之间的相互定性关系。具有以下优点：轴向功能梯度Euler-Bernoulli变截面梁屈曲临界荷载的控制方程是一组复杂的变系数常微分特征方程，研究避免了用迭代方法计算超越方程的困难和繁杂，采用变步长非均匀网格布点方式求解的梁屈曲临界荷载和已有文献计算结果完全吻合，而等步长均匀网格布点方式的计算值与实际值误差较大或甚至失真。当梁的轴向荷载由负值向正值变化，即轴向荷载由拉力转化为压力，梁的第1阶固有频率逐渐减小，当其达到屈曲临界荷载时，梁的第1阶固有频率为零。相对边界条件C-H和C-C来说，边界条件为H-H时约束最弱，其屈曲临界荷载最小，轴向荷载对第1阶固有频率的影响也越明显。