功能梯度梁靜力與動力分析的格林函數(shù)法

鄧先琪，蘇 成,2，馬海濤

(1. 華南理工大學土木與交通學院，廣州 510640；2. 華南理工大學亞熱帶建筑科學國家重點實驗室，廣州 510640；3. 廣州大學工程抗震研究中心，廣州 510405)

功能梯度材料(Functionally graded material)[1]指由兩種或兩種以上不同材料復合而成，且各組成成分的結(jié)構(gòu)和性能在空間上沿著某一方向呈連續(xù)梯度變化的一種新型非勻質(zhì)復合材料。由于其新穎的材料可設計性思想和優(yōu)良的結(jié)構(gòu)性能，近年來引起學者的廣泛關(guān)注和研究[2?5]。文獻[6 ? 8]對功能梯度材料的力學性能與研究進展進行了綜述。

目前，關(guān)于材料性能沿梁高度方向梯度變化的功能梯度梁靜動態(tài)響應分析問題，解析方法和數(shù)值方法已有大量的研究。Sankar[9]假定材料彈性模量沿梁高度方向按指數(shù)函數(shù)變化，采用Euler-Bernoulli 梁理論給出了功能梯度簡支梁在正弦分布荷載作用下的彈性力學解；而Zhong 和Yu[10]則基于Airy 應力函數(shù)給出了功能梯度懸臂梁的二維彈性力學解。Li[11]提出了功能梯度Timoshenko 梁和Euler-Bernoulli 梁的靜力與動力分析統(tǒng)一列式。Simsek[12]利用不同的高階梁理論分析了功能梯度梁的基本頻率。Chakraborty 等[13]基于一階剪切變形理論提出了一種新的梁單元來研究功能梯度梁的靜力、自由振動和波動力學行為；而Kadoli 等[14]則基于高階剪切變形理論，提出了高階梁單元來進行功能梯度梁的靜力分析。Filippi 等[15]基于一類統(tǒng)一列式提出了一種功能梯度梁單元模型，并對功能梯度梁靜力問題做了相關(guān)的研究。Pradhan和Murmu[16]采用微分求積法研究了材料參數(shù)沿高度方向呈冪函數(shù)形式分布的彈性地基功能梯度梁和層合梁的振動問題。Jiao 等[17]基于一階剪切變形與經(jīng)典梁理論，采用微分求積法分析了功能梯度梁的線性彎曲問題。李世榮和劉平[18]基于Euler-Bernoulli 梁理論，比較了功能梯度梁和均勻梁的控制微分方程，得到了功能梯度梁與均勻梁靜動態(tài)解之間的相似轉(zhuǎn)換關(guān)系。王瑄和李世榮[19]基于Levinson 高階剪切變形理論，研究了兩端簡支功能梯度梁自由振動頻率與經(jīng)典梁理論下參考勻質(zhì)梁自由振動頻率之間的解析轉(zhuǎn)換關(guān)系。以上文獻中僅考慮了材料沿梁高度梯度變化的情況。

相比于材料性能沿梁高度方向梯度變化的功能梯度梁，關(guān)于軸向功能梯度梁靜力與動力分析的研究相對較少。軸向功能梯度梁的控制微分方程本質(zhì)上是一個變系數(shù)常微分方程，通常情況下很難獲得該類方程的解析解。在數(shù)值方法方面，Murín 和Kuti?[20]采用傳遞系數(shù)描述功能梯度梁在空間上變化的彈性模量，建立了功能梯度Euler-Bernoulli 梁的等效剛度矩陣。Alshorbagy 等[21]提出了材料沿軸向或高度方向梯度變化的功能梯度梁動力特性分析的有限元法。汪亞運等[22]基于切比雪夫多項式展開技術(shù)研究了軸向功能梯度變截面梁的自由振動問題。

本文主要研究彈性模量沿軸向梯度變化的功能梯度Euler-Bernoulli 等截面梁的靜力與動力分析問題。由于這種梁的彈性模量沿軸向呈函數(shù)變化，因此采用傳統(tǒng)有限元法進行結(jié)構(gòu)分析時，不可避免地要進行有限單元網(wǎng)格細分，隨之帶來的問題是計算規(guī)模的增大與計算效率的降低。為了提高功能梯度梁靜動力問題的計算精度與效率，本文首先將功能梯度梁靜力分析的控制微分方程轉(zhuǎn)化為與勻質(zhì)材料等截面梁靜力分析控制微分方程相一致的形式，然后利用無限長勻質(zhì)材料等截面梁靜力問題的格林函數(shù)[23]，提出了功能梯度梁靜力分析的格林函數(shù)法。在此基礎上，進一步推導獲得功能梯度梁的柔度矩陣，據(jù)此建立功能梯度梁的運動方程，開展功能梯度梁的動力特性分析和動力響應分析。

1 功能梯度梁靜力分析

1.1 基本方程

考慮彈性模量沿軸向梯度變化的功能梯度梁靜力分析的控制微分方程為：

式中：x為梁的軸向坐標；u(x)和v(x)分別為梁的軸向位移和橫向位移；E為參考勻質(zhì)梁的彈性模量； η(x) 為反映彈性模量沿軸向變化的函數(shù)；A和I分別為梁的截面面積和截面慣性矩；p(x)和q(x)分別為梁的軸向和橫向分布荷載。

根據(jù)微分運算關(guān)系，式(1)可以轉(zhuǎn)化為：式中，F(xiàn)1(x)和F2(x)分別為梁的等效分布荷載，表示為：

式中：

由式(2)可以看出，功能梯度梁靜力分析的控制微分方程可以寫成與勻質(zhì)梁靜力分析的控制微分方程相似的形式，唯一的差別在于式(2)中的等效分布荷載與方程左端梁的位移相關(guān)聯(lián)。盡管如此，仍可以利用勻質(zhì)梁靜力問題的格林函數(shù)[23]進行求解。

利用勻質(zhì)梁靜力問題的格林函數(shù)進行求解時，在內(nèi)力計算上存在一些不同。考慮到材料的非勻質(zhì)性，功能梯度梁內(nèi)力與位移之間的微分關(guān)系為：

式中，b(x)見式(4)。

由式(5)可見，在功能梯度梁的靜力分析中，內(nèi)力與位移微分關(guān)系在形式上與勻質(zhì)梁問題的內(nèi)力與位移微分關(guān)系并不完全一致。因此，在采用勻質(zhì)梁靜力問題的格林函數(shù)進行求解時，需要考慮函數(shù) η(x) 及修正項b(x)M(x)對內(nèi)力的影響。

利用式(5)，可將式(3)給出的等效荷載表示為：

由式(6)可見，等效荷載可以寫成與內(nèi)力N(x)、Q(x)和M(x)相關(guān)的形式，因此在采用格林函數(shù)法分析時，需要建立關(guān)于內(nèi)力響應的補充方程，方可求解。

1.2 靜力響應分析算法

利用功能梯度梁靜力分析控制微分方程和勻質(zhì)梁靜力分析控制微分方程在形式上的相似性，可以將功能梯度梁轉(zhuǎn)化為如圖1 所示的等效勻質(zhì)梁ab 進行求解。其中，x?y為梁的坐標系，梁端點軸向坐標為xa和xb，彈性模量取為參考勻質(zhì)梁的彈性模量E，截面面積和截面慣性矩分別取為功能梯度梁的截面面積A和截面慣性矩I，軸向和橫向分布荷載分別取為式(6)所示的等效分布荷載F1(x)和F2(x)。

圖1 與功能梯度梁等效的勻質(zhì)梁Fig. 1 A uniform beam equivalently to the functionally graded beam

采用域外虛荷載方法[24]，將等效勻質(zhì)梁ab 嵌入到具有相同材料性質(zhì)與截面性質(zhì)的無限長梁中，并在梁ab 域外 ξj(j=1,2)處布設6 個大小未知的集中虛荷載Xij(i=1,2,3;j=1,2)。根據(jù)疊加原理，在域內(nèi)等效荷載Fi(x)(i=1,2)與域外虛荷載Xij(i=1,2,3;j=1, 2)的共同作用下，梁ab 上任意一點的位移與內(nèi)力響應為

理論上，此式與式(10)在同一點處應給出相同的內(nèi)力值。基于此，在梁域內(nèi)選擇(l+1)個點，利用軸力N(x) 相等可以得到(l+1)個方程；類似的，在梁域內(nèi)選擇(m+1 )個點，利用彎矩M(x)相等可以得到(m+1)個方程。最終可得：

式中，H的各行取決于式(13)系數(shù)矩陣在所選擇點處的取值，而A? 、B? 和C?的各行取決于式(10)中對應影響矩陣在所選擇點處的取值。

聯(lián)立式(12)和式(14)，可以解得：

式中：

式中，I為 6×6的單位矩陣。

從以上推導可見，由式(15)求得的基本未知向量X和Z分別包含6 個虛荷載和(l+m+2)個內(nèi)力系數(shù)，基本未知量數(shù)目與梁的子段劃分無關(guān)，即與子段數(shù)n無關(guān)。

將X和Z代入式(10)，即可得到功能梯度梁內(nèi)任意一點的位移和內(nèi)力響應為：

式中：

2 功能梯度梁動力分析

2.1 柔度矩陣

令R0=0，由式(17)可以得到由子段荷載向量f引起的梁各子段中點xk(k=1,2,···,n)處的位移為：

式中，w=[u(x1)v(x1)u(x2)v(x2) ···u(xn)v(xn)]T，f見 式(11)， δ 可 由 δ2(x) 前 兩 行 元 素 在xk(k=1, 2,···,n) 即 為 功能 梯 度 梁 對應 于w和f的 柔 度矩陣。

2.2 運動方程

阻尼矩陣可采用Rayleigh 阻尼模型生成。

2.3 動力特性和動力響應分析

式(21)可以退化為功能梯度梁無阻尼自由振動方程。利用傳統(tǒng)的特征值問題求解方法，可以獲得功能梯度梁的自振頻率和振型。

在動力荷載f(t)的作用下，采用傳統(tǒng)的振型分解法或直接積分法求解式(21)，可以獲得功能梯度梁的位移響應w(t) 、速度響應w˙(t)和加速度響應w¨(t)。同時，考慮外荷載、慣性力和阻尼力的作用，由式(17)可以進一步獲得功能梯度梁內(nèi)任意一點處的位移和內(nèi)力響應為：

3 數(shù)值算例

本文所采用的格林函數(shù)法和有限元法計算程序均采用Python 編程語言進行編寫，并在桌面電腦(Intel?Core?i7-2620 CPU, 8-GB RAM)上運行。

3.1 靜力響應分析

圖2 一端固支、一端簡支功能梯度梁Fi g. 2 A clamped-simply supported functionally graded beam

分別采用解析法、格林函數(shù)法和有限元法計算功能梯度梁的位移和內(nèi)力。根據(jù)材料力學理論，該功能梯度梁的軸向與橫向位移、軸力與彎矩解析解表達式分別為：

在格林函數(shù)法中，采用分段矩形積分公式計算式(7)中的域內(nèi)積分項，分別取n=4 、 10 、 20、50 和 100以考慮子段數(shù)對計算精度的影響。在軸向和橫向三角形分布荷載作用下，梁的軸力和彎矩分別為關(guān)于x的2 次和3 次多項式，因此分別將式(8)中的軸力和彎矩表示為相應次數(shù)的多項式，即l=2，m=3 ，并在梁內(nèi)均勻選取(l+1 )和(m+1)個點分別建立關(guān)于軸力和彎矩的補充方程。將兩端虛荷載施加位置與梁左右端點的距離表示為d1和d2，分別考慮d1=d2=0.2 m、d1=d2=0.4 m和d1=0.2 m及d2=0.4 m三種情況對格林函數(shù)法計算結(jié)果的影響。

在有限元法中，將功能梯度梁離散成Ne個基于Hermite 插值的平面二節(jié)點梁單元，分別取Ne=4 、 10 、 20 、 50 和 100以 考 慮 離 散 單 元 數(shù) 對計算精度的影響，其中每個梁單元的彈性模量取為梁單元中點處的彈性模量值。

表1 給出了功能梯度梁跨中點C 的軸向位移和橫向位移結(jié)果，其收斂情況如圖3 所示。表2給出了功能梯度梁固定端點A 的軸力和彎矩結(jié)果，其收斂情況如圖4 所示。由上述圖表可見，當子段數(shù)n≥20時，格林函數(shù)法的計算結(jié)果與解析解高度吻合，說明所提出方法具有理想的計算精度。此外，由計算結(jié)果還可以看出，虛荷載位置的三種選擇對格林函數(shù)法計算結(jié)果沒有任何影響，說明所提出方法具有良好的數(shù)值穩(wěn)定性。

表1 功能梯度梁跨中點C 的軸向與橫向位移Table 1 Axial and transverse displacement of node C at midspan of the functionally graded beam

圖3 功能梯度梁跨中點C 的軸向與橫向位移Fig. 3 Axial and transverse displacement of node C at midspan of the functionally graded beam

表2 功能梯度梁固定端點A 的軸力與彎矩Table 2 Axial force and bending moment of node A at fixed end of the functionally graded beam

圖4 功能梯度梁固定端點A 的軸力和彎矩Fig. 4 Axial force and bending moment of node A at fixed end of the functionally graded beam

在計算效率方面，格林函數(shù)法中的基本未知量為虛荷載X和內(nèi)力系數(shù)Z，詳見式(12)和式(14)，未知量數(shù)目并不受子段數(shù)n的影響，因此即使n取較大值也不會影響方程組的求解效率。而在有限元法中，隨著離散單元數(shù)Ne的增加，基本未知量數(shù)目和剛度矩陣規(guī)模也相應增多，從而計算工作量也會顯著增加，計算時間延長。

3.2 動力特性分析

考慮圖2 所示一端固支、另一端簡支的功能梯度梁。梁的截面面積、截面慣性矩和長度與算例3.1 中相應參數(shù)相同。假設梁的彈性模量和密度均沿軸向呈3 次多項式函數(shù)形式變化，表達式為：

式 中：EL=210 GPa和ER=390 GPa分 別 表 示 梁左、右兩端處彈性模量的取值；ρL=7800 kg/m3和ρR=3960 kg/m3分別表示梁左、右兩端處密度的取值。

分別采用格林函數(shù)法和有限元法計算功能梯度梁的自振頻率。在格林函數(shù)法中，反映虛荷載位置的距離取為d1=d2=0.2 m。為了考察子段數(shù)n對計算精度的影響，分別取n=10 、 20 、 30、40 、 50 和 200。由于梁在自由振動過程中其內(nèi)力分布不再滿足簡單的多項式分布形式，為了考察彎矩多項式次數(shù)m對橫向自由振動固有頻率的影響，分別取m=4、6、8 和10，而軸力多項式次數(shù)則取為l=10。在有限元法中，將功能梯度梁離散成Ne個基于Hermite 插值的平面二節(jié)點梁單元，采用一致質(zhì)量矩陣。分別取Ne=10 、 20、30 、 40 、 50 和 200 以考察離散單元數(shù)Ne對計算精度的影響，其中每個梁單元的彈性模量取為梁單元中點處的彈性模量值。

表3 和圖5 給出了功能梯度梁的前三階橫向自由振動固有頻率的計算結(jié)果。由計算結(jié)果可見，當m>6 時，在子段數(shù)n和離散單元數(shù)Ne取為相同的情況下，格林函數(shù)法計算得到的前三階橫向自由振動固有頻率總是比有限元法計算得到的結(jié)果精度更高，收斂速度更快。也就是說，當彎矩多項式次數(shù)足以反映真實彎矩分布時，要達到相同的精度，格林函數(shù)法所需要的子段數(shù)n少于有限元法所需要的離散單元數(shù)Ne，因此其計算效率更高。

表3 功能梯度梁前三階橫向自由振動固有頻率Table 3 The first three natural frequencies of functionally graded beam’s transverse free vibration

圖5 功能梯度梁前三階橫向自由振動固有頻率Fig. 5 The first three natural frequencies of functionally graded beam’s transverse free vibration

3.3 動力響應分析

考慮圖2 所示一端固支、另一端簡支的功能梯度梁，其幾何參數(shù)和材料參數(shù)與算例3.2 中的相應參數(shù)相同。梁受橫向均布突加荷載q(t)作用，荷載隨時間變化曲線如圖6 所示。選取Rayleigh 阻尼模型，阻尼比取為 ζ=0.05(對應的自振頻率ω1=1324.42 rad/s和ω2=4317.06 rad/s)。

圖6 荷載-時間曲線Fig. 6 Load-time curve

分別采用格林函數(shù)法和有限元法計算功能梯度梁的位移、速度和加速度。在格林函數(shù)法中，反映虛荷載位置的距離取為d1=d2=0.2 m，內(nèi)力多項式次數(shù)取為l=m=8 。為了考察子段數(shù)n對計算精度的影響，分別取n=10 、 20 和 30。在有限元法中，將功能梯度梁離散成Ne個基于Hermite插值的平面二節(jié)點梁單元，采用一致質(zhì)量矩陣。取Ne=200以獲得相對精確解，其中每個梁單元的彈性模量取為梁單元中點處的彈性模量值。兩種方法均采用Newmark- β積分格式求解運動方程，積分步長均取為 ?t=1.0×10?4s ，共計算 500個時間步。

圖7 給出了功能梯度梁跨中點C 的橫向位移、速度和加速度結(jié)果，上述結(jié)果在0.03 s 后的時間內(nèi)基本不再變化。從圖中可見，當子段數(shù)n=30時，格林函數(shù)法計算得到的跨中點C 的橫向位移、速度和加速度響應與有限元法結(jié)果高度吻合，再次說明所提出方法具有理想的計算精度。

圖7 功能梯度梁跨中點C 橫向位移、速度和加速度Fig. 7 Transverse displacement, velocity, acceleration of node C at midspan of the functionally graded beam

4 結(jié)論

本文提出了一種軸向功能梯度材料Euler-Bernoulli 梁靜力與動力分析的格林函數(shù)法。將功能梯度梁靜力分析的控制微分方程轉(zhuǎn)化為與勻質(zhì)材料梁靜力分析控制微分方程相一致的形式，從而可以利用勻質(zhì)材料梁靜力問題的格林函數(shù)進行求解。在此基礎上，進一步可以推導得到功能梯度梁的柔度矩陣，并建立功能梯度梁的運動方程，開展功能梯度梁動力特性分析與動力響應分析。數(shù)值算例表明，對于功能梯度梁的靜力和動力問題，格林函數(shù)法比有限元法具有更高的計算精度和計算效率。