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

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

(安徽工程大學 建筑工程學院，安徽 蕪湖 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階固有頻率的影響也越明顯。