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

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

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

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

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

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

1 功能梯度梁靜力分析

1.1 基本方程

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

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

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

式中：

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

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

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

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

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

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

1.2 靜力響應(yīng)分析算法

利用功能梯度梁靜力分析控制微分方程和勻質(zhì)梁靜力分析控制微分方程在形式上的相似性，可以將功能梯度梁轉(zhuǎn)化為如圖1 所示的等效勻質(zhì)梁ab 進(jìn)行求解。其中，x?y為梁的坐標(biāo)系，梁端點(diǎn)軸向坐標(biāo)為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)處布設(shè)6 個(gè)大小未知的集中虛荷載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 上任意一點(diǎn)的位移與內(nèi)力響應(yīng)為

理論上，此式與式(10)在同一點(diǎn)處應(yīng)給出相同的內(nèi)力值?；诖耍诹河騼?nèi)選擇(l+1)個(gè)點(diǎn)，利用軸力N(x) 相等可以得到(l+1)個(gè)方程；類似的，在梁域內(nèi)選擇(m+1 )個(gè)點(diǎn)，利用彎矩M(x)相等可以得到(m+1)個(gè)方程。最終可得：

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

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

式中：

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

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

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

式中：

2 功能梯度梁動(dòng)力分析

2.1 柔度矩陣

令R0=0，由式(17)可以得到由子段荷載向量f引起的梁各子段中點(diǎn)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) 即 為 功能 梯 度 梁 對(duì)應(yīng) 于w和f的 柔 度矩陣。

2.2 運(yùn)動(dòng)方程

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

2.3 動(dòng)力特性和動(dòng)力響應(yīng)分析

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

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

3 數(shù)值算例

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

3.1 靜力響應(yīng)分析

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

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

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

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

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

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

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

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

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

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

3.2 動(dòng)力特性分析

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

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

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

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

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

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

3.3 動(dòng)力響應(yīng)分析

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

圖6 荷載-時(shí)間曲線Fig. 6 Load-time curve

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

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

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

4 結(jié)論

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