功能梯度梁臨界荷載與固有頻率的微分求積法計算

熊海超，葛仁余，馬國強，劉小雙，余本源

(安徽工程大學(xué) 建筑工程學(xué)院，安徽 蕪湖 241000)

功能梯度材料是一種新型材料，它是將多種性能各異的材料按照設(shè)計意愿形成性能連續(xù)變化的合成材料，以滿足各種特殊工程結(jié)構(gòu)的需要。但是，由于功能梯度材料的物理性質(zhì)連續(xù)變化，給其力學(xué)分析帶來一定的難題，除非一些特殊情況，一般用解析方法分析功能梯度材料的振動問題十分困難，如文獻(xiàn)[1]用能量法獲得了錐形變截面樁屈曲臨界荷載的解析解。目前，數(shù)值方法是用來研究軸向功能梯度變截面梁的自由振動和屈曲臨界荷載的主要手段，文獻(xiàn)[2-3]運用精確的動態(tài)剛度方法，獲得了Euler-Bernoulli梁和Timoshenko梁的固有頻率和屈曲臨界荷載；文獻(xiàn)[4]運用微分求積單元法研究了軸向功能梯度Euler-Bernoulli梁自由振動和屈曲臨界荷載的問題；基于Euler-Bernoulli梁理論和Timoshenko梁理論，文獻(xiàn)[5-6]運用插值矩陣法研究了軸向功能梯度變截面梁自由振動和屈曲臨界荷載問題。將邊界條件轉(zhuǎn)換為簡便的形式，避免了剛度系數(shù)矩陣的無窮大值，文獻(xiàn)[7]研究了軸向功能梯度變截面彈性地基梁的屈曲臨界荷載問題。文獻(xiàn)[8]采用局部微分求積法研究了帶有彈性約束的軸向功能梯度變截面樁屈曲臨界荷載問題。利用均勻梁單元的形狀函數(shù)，提出了一種新的梁單元。文獻(xiàn)[9]采用有限元計算方法，研究了軸向功能梯度梁的自由振動和穩(wěn)定性問題。

研究采用切比雪夫多項式的根作為非均勻網(wǎng)格布點方式，由微分求積法(DQM)研究受軸向荷載作用的功能梯度材料Euler-Bernoulli變截面梁的屈曲臨界荷載、固有頻率，以及軸向荷載對固有頻率的影響。基于Euler-Bernoulli梁理論，獲得求解功能梯度材料Euler-Bernoulli梁的屈曲臨界荷載和固有頻率的變系數(shù)常微分方程，再根據(jù)微分求積法理論，將變系數(shù)常微分方程的特征值問題轉(zhuǎn)化為線性代數(shù)方程組的特征值問題，可一次性地獲得梁的屈曲臨界荷載、固有頻率，以及軸向荷載對梁的第1階固有頻率的影響。

1 計算模型

1.1 Euler-Bernoulli梁穩(wěn)定性的基本理論

圖1 軸向功能梯度變截面Euler-Bernoulli梁

軸向功能梯度變截面Euler-Bernoulli梁如圖1所示。功能梯度材料Euler-Bernoulli變截面梁的長度為

l

，材料性能和截面面積沿軸向

x

任意連續(xù)變化。設(shè)梁自由振動時的撓度為

w(x,t)

，材料的彈性模量為

E(x)=E

f

(x)

，材料的密度為

ρ(x)=ρ

f

(x)

，截面面積為

A(x)=A

f

(x)

，截面轉(zhuǎn)動慣量為

I(x)=I

f

(x)

，它們均為

x

的連續(xù)函數(shù)。其中，

E

、

ρ

、

A

、

I

對應(yīng)于變截面梁在左端邊界

x=

0位置處材料的彈性模量、密度、截面面積和截面慣性矩；

f

(x)

為

x

的函數(shù)

(i=

1

,

2

,

3

,

4

)

；設(shè)軸向荷載

P>

0為壓力，

P<

0為拉力，則在軸向荷載

P

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

(

1

)

研究主要考慮梁的諧波振動，則有：

w(x,t)=W(x)

sin

ωt,

(

2

)

式中，

ω

為梁振動角頻率

;t

為時間變量

;w(x,t)

為梁撓度

;W(x)

為梁振型函數(shù)。聯(lián)立式

(

1

)

和式

(

2

)

，可得軸向功能梯度Euler-Bernoulli變截面梁在軸向荷載

P

作用下的振動控制方程：

(

3

)

(

4

)

將式

(

4

)

代入式

(

3

)

，可得軸向荷載

P

作用下軸向功能梯度Euler-Bernoulli變截面梁的自由振動控制方程：

(

5

)

當(dāng)軸向荷載

P

達(dá)到屈曲臨界荷載時，功能梯度材料Euler-Bernoulli變截面梁的固有頻率

ω=

0，將式

(

5

)

轉(zhuǎn)化為式

(

6

)

關(guān)于屈曲臨界荷載

λ

的變系數(shù)常微分方程的特征值問題：

(

6

)

1.2 Euler-Bernoulli梁的邊界條件

Euler-Bernoulli梁的轉(zhuǎn)角

θ

，彎矩

M

和剪力

T

分別為：

(

7

)

簡支-簡支梁

(

H-H

)

的邊界條件是兩端撓度和彎矩為零，即

W(

0

)=

0

,W

(

0

)=

0

;W(

1

)=

0

,W

(

1

)=

0

,

(8a)

固端-固端梁

(

C-C

)

的邊界條件是兩端撓度和轉(zhuǎn)角為零，即

W(

0

)=

0

,W

(

0

)=

0

;W(

1

)=

0

,W

(

1

)=

0

,

(8b)

固端-簡支梁

(

C-H

)

的邊界條件是固端撓度和轉(zhuǎn)角為零，簡支端撓度和彎矩為零，即

W(

0

)=

0

,W

(

0

)=

0

;W(

1

)=

0

,W

(

1

)=

0

,

(8c)

因此，功能梯度材料Euler-Bernoulli變截面梁固有頻率的計算實際成為求解邊界條件式

(

8

)

下變系數(shù)常微分方程式

(

5

)

的特征值問題；屈曲臨界荷載的計算實際成為求解邊界條件式

(

8

)

下變系數(shù)常微分方程式

(

6

)

的特征值問題。研究采用微分求積法對上述問題展開數(shù)值方法分析。

2 微分求積法數(shù)值計算方法

功能梯度材料Euler-Bernoulli變截面梁微分求積法計算模型如圖2所示。研究采用切比雪夫多項式的根作為節(jié)點坐標(biāo)的方式如式(9a)所示，均勻等步長分布節(jié)點坐標(biāo)的方式如式(9b)所示：

圖2 軸向功能梯度材料Euler-Bernoulli梁的網(wǎng)格劃分模型

(

9a

)

(

9b

)

考慮一維函數(shù)

W(ξ)

在區(qū)間

[

0

,

1

]

上可微，將區(qū)間

[

0

,

1

]

劃分為

n

段，0

=ξ

，

ξ

，

ξ

，…

,ξ

-1，

ξ

=

1，在每個節(jié)點上函數(shù)及其導(dǎo)數(shù)用區(qū)間

[

0

,

1

]

上

n

+1個節(jié)點函數(shù)值的加權(quán)線性和近似表示。這里將梁的振型函數(shù)

W(ξ)

用拉格朗日

(

Lagrange

)

插值函數(shù)來描述，即

(

10

)

式中，

l

(ξ)

為拉格朗日多項式，其表達(dá)式為：

(

11

)

由式

(

10

)

對函數(shù)

W(ξ)

求一階導(dǎo)數(shù)，得到：

(

12

)

將式

(

12

)

在

n

段區(qū)間

[

0

,

1

]

上進(jìn)行離散，從而進(jìn)一步得到：

(

13

)

依次類推可得：

(

14

)

由式

(

14

)

可知：

(

15

)

(

16

)

(

17

)

由于

(

18

)

因此，各階導(dǎo)數(shù)的加權(quán)系數(shù)矩陣之間的關(guān)系為：

(

19

)

將式

(

6

)

計算屈曲臨界荷載的微分方程組中變系數(shù)寫成對角陣向量形式，

I

為單位陣，即

(

20

)

將式

(

6

)

微分方程組寫成代數(shù)方程組向量形式為：

C

W

-

λC

W=

0

,

(

21

)

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

(

C-C

)

情況為例進(jìn)行討論，則相應(yīng)的邊界條件向量形式可表示為：

(

22

)

這里，

I

是0和

n

，

[

…

]

為矩陣

[

…

]

的第

I

行元素。研究求解如式

(

6

)

的四階微分方程時，運用節(jié)點替代法

(δ

法

)

來處理邊界條件，即將邊界條件利用與梁兩端相隔

δ

的兩個點來替代

(

10

<δ<

10

)

，取節(jié)點序列為

{ξ}={

0

,δ,

…

,ξ

,

…

,

1-

δ,

1

}

(j=

3

,

4

,

…

,n

-2

),

(

23

)

(

24

)

3 數(shù)值算例和問題討論

設(shè)軸向功能梯度Euler-Bernoulli變截面梁的長度為

l

，截面形狀為矩形，材料的彈性模量和密度分別為：

梁的截面面積和慣性矩分別為：

式中，

c

和

c

分別是梁截面的寬、高錐度系數(shù)

;x

是以梁的左端為起點沿軸線方向的橫坐標(biāo)；當(dāng)

c

=c

=

0時為等截面梁。運用微分求積法計算軸向功能梯度Euler-Bernoulli變截面梁屈曲臨界荷載的變系數(shù)常微分方程式，即式

(

6

)

，離散單元數(shù)取

n=

20，非均勻網(wǎng)格

(

切比雪夫多項式根作為節(jié)點

)

和均勻網(wǎng)格兩種布點方式的微分求積法計算值對比如表1所示。由表1計算結(jié)果可知，切比雪夫多項式根作為布點方式的計算值與文獻(xiàn)

[

9

]

結(jié)果吻合良好，而均勻等步長布點方式的計算值與文獻(xiàn)

[

9

]

結(jié)果誤差較大甚至失真。大量的數(shù)值實驗研究表明，微分求積法求解變系數(shù)微分方程時，采用式(9

b

)均勻節(jié)點劃分網(wǎng)格時，微分求積法計算數(shù)值不穩(wěn)定甚至失真，而一些非均勻節(jié)點劃分網(wǎng)格往往能得到很好的結(jié)果，如切比雪夫多項式的根作為離散節(jié)點分布形式(式(9

a

))。

表1 在不同布點方式下，軸向功能梯度變截面梁的無量綱屈曲臨界荷載計算值比對

在C-C、C-H與H-H不同邊界條件和梁截面錐度系數(shù)下，以切比雪夫多項式根作為變步長非均勻網(wǎng)格布點方式，取

n=

20，微分求積法求解變系數(shù)常微分方程式

(

5

)

獲得的荷載

λ

與第1階固有頻率

μ

的關(guān)系曲線圖如圖3所示。由圖3可知：①隨著軸向荷載由負(fù)值向正值增大，即軸向荷載由拉力向壓力轉(zhuǎn)化，梁的第1階固有頻率逐漸減小；當(dāng)軸向荷載較小時，梁振動頻率近似呈線性降低，當(dāng)軸向荷載即將接近屈曲臨界荷載值時，梁的第1階橫向振動頻率開始非線性急劇降低；②當(dāng)邊界條件為H-H時，梁的邊界約束最弱，軸向荷載對第1階固有頻率的影響明顯，屈曲臨界荷載最小；當(dāng)邊界條件為C-C時，梁的邊界約束最強，屈曲臨界荷載最大；③隨著壓力的增加，梁第1階固有頻率逐漸趨于零，當(dāng)?shù)?階固有頻率為零時，對應(yīng)的軸向荷載即為梁的屈曲臨界荷載，在邊界條件C-C、C-H和H-H下，截面寬錐度系數(shù)

c

=

0.2和截面高錐度系數(shù)

c

=

0.2時，梁的無量綱屈曲臨界荷載

λ

分別為37.602 327 3、19.335 410 3、9.597 120 84；截面寬錐度系數(shù)

c

=

0.8和截面高錐度系數(shù)

c

=

0.8時，無量綱屈曲臨界荷載

λ

分別為21.781 289 8、11.218 322 3、5.622 775 59；截面寬錐度系數(shù)

c

=0.4和截面高錐度系數(shù)

c

=

0.4時，無量綱屈曲臨界荷載

λ

分別為0.664 879 02、1.371 282 41、0.707 535 808，微分求積法求解變系數(shù)常微分方程式

(

5

)

的計算值與表1的計算值自行吻合。

表2 軸向功能梯度材料Euler-Bernoulli變截面梁的無量綱屈曲臨界荷載計算值

圖3 軸向荷載對Euler-Bernoulli梁第1階固有頻率的影響示意圖

4 結(jié)論

研究基于Euler-Bernoulli梁理論，建立了求解功能梯度材料變截面梁的屈曲臨界荷載和固有頻率的變系數(shù)常微分方程，再運用微分求積法理論，通過縝密的數(shù)學(xué)推導(dǎo)，將變系數(shù)常微分方程的特征值問題轉(zhuǎn)化為線性代數(shù)方程組的特征值問題，一次性地獲得梁的屈曲臨界荷載以及軸向力與梁固有頻率之間的相互定性關(guān)系。具有以下優(yōu)點：軸向功能梯度Euler-Bernoulli變截面梁屈曲臨界荷載的控制方程是一組復(fù)雜的變系數(shù)常微分特征方程，研究避免了用迭代方法計算超越方程的困難和繁雜，采用變步長非均勻網(wǎng)格布點方式求解的梁屈曲臨界荷載和已有文獻(xiàn)計算結(jié)果完全吻合，而等步長均勻網(wǎng)格布點方式的計算值與實際值誤差較大或甚至失真。當(dāng)梁的軸向荷載由負(fù)值向正值變化，即軸向荷載由拉力轉(zhuǎn)化為壓力，梁的第1階固有頻率逐漸減小，當(dāng)其達(dá)到屈曲臨界荷載時，梁的第1階固有頻率為零。相對邊界條件C-H和C-C來說，邊界條件為H-H時約束最弱，其屈曲臨界荷載最小，軸向荷載對第1階固有頻率的影響也越明顯。