徐 超，王 騰

(西北工業(yè)大學(xué) 航天學(xué)院，西安 710072)

功能梯度材料的宏觀材料特性表現(xiàn)出空間連續(xù)變化的性質(zhì)。材料成分和性質(zhì)的變化會(huì)改變彈性波在結(jié)構(gòu)中的傳播行為，研究彈性波在功能梯度等非均勻材料結(jié)構(gòu)中的傳播特性對(duì)發(fā)展新材料結(jié)構(gòu)的損傷檢測(cè)和健康監(jiān)測(cè)方法具有重要的意義［1］。

與傳統(tǒng)均質(zhì)材料不同，功能梯度材料的彈性模量、泊松比等宏觀性質(zhì)都隨空間位置而變化，采用解析解或半解析方法研究功能梯度材料結(jié)構(gòu)中的波傳播行為十分困難。文獻(xiàn)［2］基于Hamilton變分原理和高階剪切變形理論推導(dǎo)了在點(diǎn)沖擊下無(wú)限大功能梯度材料板應(yīng)力波的近似解析解。文獻(xiàn)［3］采用一階剪切變形理論和小應(yīng)變的應(yīng)變－位移關(guān)系，研究了四邊固支板中，不同功能梯度材料指數(shù)對(duì)波傳播行為的影響。文獻(xiàn)［4］采用連續(xù)材料模型，利用Legendre多項(xiàng)式展開方法獲得了功能梯度材料板的頻散和功率流解。文獻(xiàn)［5］應(yīng)用冪級(jí)數(shù)展開方法研究了Lamb波在功能梯度材料板結(jié)構(gòu)中的傳播行為。然而，這些工作多局限于一維梯度材料，并且研究對(duì)象均為邊界條件和幾何較為簡(jiǎn)單的結(jié)構(gòu)。

譜單元法是求解復(fù)雜結(jié)構(gòu)彈性波傳播行為的一種新型數(shù)值方法。文獻(xiàn)［6］最早在計(jì)算流體問(wèn)題研究中提出了譜單元法，其充分結(jié)合了譜方法的快速收斂和有限元法對(duì)復(fù)雜結(jié)構(gòu)適應(yīng)性好的優(yōu)點(diǎn)，特別適合求解瞬變和高頻動(dòng)力學(xué)問(wèn)題。文獻(xiàn)［7］將譜單元推廣用于彈性波傳播的模擬，并求解了一維均勻材料桿結(jié)構(gòu)中的波傳播問(wèn)題。文獻(xiàn)［8－10］則進(jìn)一步地將譜單元法用于二維和三維結(jié)構(gòu)波傳播問(wèn)題。由于功能梯度材料性質(zhì)在結(jié)構(gòu)中的空間變化特性，應(yīng)用譜單元法求解功能梯度材料的結(jié)構(gòu)波傳播問(wèn)題的研究還不多。文獻(xiàn)［11－12］采用頻域譜單元建立了新的功能梯度材料梁?jiǎn)卧?，用于沖擊載荷下功能梯度梁的波傳播問(wèn)題。文獻(xiàn)［12］建立了一種基于切比雪夫多項(xiàng)式時(shí)域譜單元，并將其用于求解平面功能梯度材料結(jié)構(gòu)中的波傳播問(wèn)題。

本文推導(dǎo)了一種具有任意四邊形形狀的Gauss－Lobatto－Legendre多項(xiàng)式時(shí)域譜單元用于功能梯度材料結(jié)構(gòu)中彈性波傳播行為的模擬。分別采用連續(xù)材料模型、分層離散模型和均勻化模型建模功能梯度材料宏觀性質(zhì)空間變化特性。將數(shù)值計(jì)算結(jié)果與理論解進(jìn)行了對(duì)比以驗(yàn)證單元有效性，研究了面內(nèi)受水平激勵(lì)的功能梯度材料板結(jié)構(gòu)的波傳播特性，對(duì)比了三種材料建模方法的差異。

1 任意四邊形譜單元推導(dǎo)

1.1 位移、形狀插值

如圖1所示，與傳統(tǒng)有限元不同，自然坐標(biāo)系下標(biāo)準(zhǔn)譜單元的節(jié)點(diǎn)非等距地分布于區(qū)域Λ∈［－1，1］2內(nèi)，單元上節(jié)點(diǎn)坐標(biāo)位置根據(jù)特殊多項(xiàng)式的零點(diǎn)確定。非等距節(jié)點(diǎn)布置能有效克服等距節(jié)點(diǎn)布置在插值時(shí)引起的龍格效應(yīng)。這里采用 Gauss－Lobatto－Legendre多項(xiàng)式確定單元中各節(jié)點(diǎn)的坐標(biāo)位置。以ξ方向?yàn)槔珿auss－Lobatto－Legendre 多項(xiàng)式為

式中:Pn－1(ξ)為n－1階勒讓德多項(xiàng)式。單元階次不同，對(duì)應(yīng)的節(jié)點(diǎn)坐標(biāo)位置不同，如當(dāng)n=5時(shí)，節(jié)點(diǎn)ξ方向坐標(biāo)為

圖1 2D譜單元及5×5節(jié)點(diǎn)分布Fig.1 2D spectral element with 5 ×5 nodes

同理，單元η方向上的節(jié)點(diǎn)坐標(biāo)也可類似確定。

仍采用拉格朗日插值函數(shù)對(duì)單元內(nèi)位移場(chǎng)進(jìn)行插值，即

式中:u，v分別為 x，y方向的位移。ψm(ξ，η)=φi(ξ)φj(η)，m=1，2，...，n2，φ 為單變量拉格朗日插值函數(shù)。

實(shí)際應(yīng)用中，剖分離散后的單元幾何形狀可能不是標(biāo)準(zhǔn)的邊長(zhǎng)為1的正方形?？紤]單元形狀為任意四邊形的情況，需進(jìn)行形狀插值建立任意四邊形單元與標(biāo)準(zhǔn)譜單元的幾何映射關(guān)系。

設(shè)任意四邊形四個(gè)角節(jié)點(diǎn)的坐標(biāo)為(xk，yk)，k=1，2，3，4，采用一階拉格朗日插值函數(shù)建立單元形狀插值關(guān)系:

1.2 功能梯度材料模型

考慮材料性質(zhì)一維變化的功能梯度材料，其彈性模量E、泊松比ν和密度ρ都沿某一方向規(guī)律變化。對(duì)于二維情況，不失一般性，假設(shè)材料性質(zhì)沿y方向變化。由于功能梯度材料的特點(diǎn)，進(jìn)行數(shù)值計(jì)算時(shí)需要采用特殊的材料模型。常用材料建模方法有三種:

(1)均勻化模型:忽略掉材料性質(zhì)的梯度變化，采用材料性能的平均值作為均勻化的材料屬性。

(2)分層離散模型:將功能梯度材料沿材料性質(zhì)變化方向離散為若干材料層，材料層內(nèi)各單元材料性質(zhì)都為同一常數(shù)。原則上講，分層模型解會(huì)隨著分層數(shù)增加趨于解析解，相應(yīng)的計(jì)算效率也會(huì)下降。

(3)連續(xù)材料模型:將材料屬性考慮為坐標(biāo)的連續(xù)函數(shù)，則單元內(nèi)材料性質(zhì)也是坐標(biāo)的連續(xù)函數(shù)，材料密度矩陣μe和彈性矩陣De可表示為:

一般地，采用體積分?jǐn)?shù)函數(shù)g(y)來(lái)描述材料屬性的變化。這里設(shè)f(y)為空間某一材料屬性值，在單元內(nèi)部，則由式(3)可知有f(y(η))。假設(shè)fu，fd分別為結(jié)構(gòu)上下表面的材料屬性值，則有

g(y)一般為冪函數(shù)形式，即

式中:yu，yd代表結(jié)構(gòu)上下界面的坐標(biāo)，h=yu－yd。χ表示材料分布的不同梯度，當(dāng)χ=0為均勻材料。

本文將分別采用上述材料建模方法研究功能梯度材料結(jié)構(gòu)中的波傳播特性，對(duì)比不同的材料模型對(duì)結(jié)構(gòu)中波傳播行為的影響。

1.3 單元?jiǎng)偠群唾|(zhì)量矩陣

利用上述位移和形狀插值函數(shù)，根據(jù)Hamilton原理建立變分方程，可推導(dǎo)獲得單元質(zhì)量矩陣Me，單元?jiǎng)偠染仃嘖e和單元上的等效外力Fe表達(dá)式為:

式中:Be矩陣為幾何矩陣，可以由位移－應(yīng)變幾何關(guān)系求出:

式(7)中 μe(η)，De(η)如為連續(xù)材料模型時(shí)采用式(4)－(6)的，其值由積分點(diǎn)位置確定;如為材料模型Ⅰ和Ⅱ時(shí)，取為相應(yīng)均勻化后的常數(shù)值。

式(7)中J為雅克比矩陣，描述任意四邊形單元與標(biāo)準(zhǔn)譜單元的幾何關(guān)系，根據(jù)式(3)有

對(duì)式(7)采用 Gauss－Lobatto－Legendre(GLL)數(shù)值積分求解，積分系數(shù)ωi＞0，可由下式確定:

由于勒讓德多項(xiàng)式的正交性，單元質(zhì)量矩陣為對(duì)角矩陣，這可顯著降低波動(dòng)力響應(yīng)時(shí)域求解的計(jì)算耗費(fèi)，較切比雪夫多項(xiàng)式時(shí)域譜單元有明顯優(yōu)勢(shì)。

2 數(shù)值算例

2.1 中心受垂直沖擊載荷的無(wú)限大FGM板

為了驗(yàn)證所推導(dǎo)單元的有效性，考慮一中心受載的無(wú)限大FGM板。文獻(xiàn)［2］采用彈性板殼理論，考慮高階剪切變形與轉(zhuǎn)動(dòng)慣量的影響，推導(dǎo)獲得了極坐標(biāo)下的解析解。本文采用平面應(yīng)變模型，利用前文推導(dǎo)的連續(xù)材料模型譜單元離散板結(jié)構(gòu)，將獲得的數(shù)值解與解析解進(jìn)行對(duì)比。

如圖2所示，板厚度為0.02 m，材料屬性沿厚度方向線性變化，即χ=1。板上、下界面材料屬性分別為鋁和氧化鋯，彈性模量分別為70GPa、151GPa，密度分別為 2700 kg/m3、3000 kg/m3，泊松比為常值 0.3。沖擊激勵(lì)形式為5個(gè)周期經(jīng)過(guò)調(diào)制的正弦信號(hào)。

圖2 板的二維模型Fig.2 Two－dimensional model of plate

應(yīng)用譜單元數(shù)值求解時(shí)，單元尺寸為0.01×0.067 m，每個(gè)單元節(jié)點(diǎn)數(shù)為5×5，積分步長(zhǎng)取2×10－8s，監(jiān)測(cè)距加載處x=2 m處的響應(yīng)和解析解對(duì)比。采用文獻(xiàn)［2］中的無(wú)量綱方法，無(wú)量綱后的時(shí)域響應(yīng)如圖3(a)、(b)所示。其中，圖3(a)給出了x方向的位移響應(yīng)，圖3(b)給出了y方向的位移響應(yīng)。由圖可知，數(shù)值解和解析解吻合較好，驗(yàn)證了本文推導(dǎo)單元的有效性。需要說(shuō)明的是，這里是依據(jù)平面應(yīng)變假設(shè)進(jìn)行數(shù)值求解，給出了近似解，目的是驗(yàn)證單元的波傳播描述能力。由于連續(xù)材料模型的單元精度也受單元規(guī)模和積分階次制約，給出的解也是數(shù)值近似解，因此圖3種仍存在小量誤差。

圖3 在x=2 m處，板上表面X方向的位移和中面Y方向的位移Fig.3 The displacements x on top surface and deflection y on middle plane of plate at x=2 m

2.2 面內(nèi)一端受水平激勵(lì)的FGM板

2.2.1 模型描述

如圖4所示，一個(gè)陶瓷－鉻混合功能梯度板在xy面內(nèi)的截面，其材料屬性見表1。板厚20mm，長(zhǎng)2 m，在左端施加水平的均勻分布載荷，激勵(lì)信號(hào)為中心頻率為50 kHz經(jīng)漢寧窗調(diào)制的正弦信號(hào)，時(shí)長(zhǎng)三個(gè)周期，信號(hào)峰值109N/m，見圖5。采用7×7節(jié)點(diǎn)譜單元計(jì)算，單元尺寸為0.005×0.05 m。動(dòng)力學(xué)方程采用中心差法求解，積分步長(zhǎng)為2×10－2μs。

圖4 板xy截面及計(jì)算網(wǎng)格Fig.4 Plane xy of plate and computational grids

圖5 中心頻率為50 kHz的激勵(lì)波形Fig.5 The excitation with the center frequency 50 kHz

公開文獻(xiàn)中對(duì)于FGM板時(shí)域求解的方法仍主要集中于分層取均值等分層離散方法［14］，本文采取三種材料模型進(jìn)行對(duì)比求解。三種材料模型采用同一網(wǎng)格，同一單元階次，同一積分方法，只是單元材料屬性模型不同，從而通過(guò)波場(chǎng)、時(shí)域響應(yīng)，相速度頻散對(duì)比，研究其對(duì)功能梯度材料結(jié)構(gòu)波傳播行為的描述能力差異。

表1 FGM材料屬性Tab.1 Functional Graded Materials properties

2.2.2 波場(chǎng)及時(shí)域響應(yīng)

左端施加激勵(lì)后，導(dǎo)波會(huì)沿著板截面?zhèn)鞑?，同時(shí)受到上下邊界的約束。

圖6給出了采用三種不同的材料模型計(jì)算出的結(jié)構(gòu)在200μs時(shí)x、y方向的位移場(chǎng)云圖。在二維結(jié)構(gòu)中，由于泊松效應(yīng)，水平方向的激勵(lì)在水平和垂直方向都引起結(jié)構(gòu)的振動(dòng)。在水平方向，能量主要集中在以對(duì)稱模式S0波傳播的模態(tài)上;在垂直方向，能量主要集中以反對(duì)稱模式A0波傳播的模態(tài)上。

由圖6可以看出，在水平和垂直方向，采用三種材料模型計(jì)算得到波場(chǎng)云圖并不相同。在x方向，200μs時(shí)，采用連續(xù)材料模型計(jì)算得到的波前位置最靠后，表明計(jì)算的波速較采用均勻模型和分散離散模型低。特別注意的是，采用均勻化材料模型只能計(jì)算得到對(duì)稱模式的波包，無(wú)法計(jì)算出緊隨其后的反對(duì)稱模式波包。在y方向，也有類似的結(jié)果。

為了進(jìn)一步定量研究不同模型對(duì)波動(dòng)響應(yīng)的影響，取板上表面的中點(diǎn)作為數(shù)值結(jié)果時(shí)域響應(yīng)的提取點(diǎn)(因一般壓電傳感器貼在結(jié)構(gòu)上表面，監(jiān)測(cè)局部響應(yīng))，給出其0～400μs的位移響應(yīng)，結(jié)果如圖7所示。

圖7(a)、(b)分別給出了監(jiān)測(cè)點(diǎn)x、y方向0～400 μs的位移響應(yīng)。由圖可知，分層離散模型和連續(xù)模型得到的縱波和橫波均含有兩個(gè)波包，即對(duì)稱和反對(duì)稱模態(tài)的波，其中縱波幅值大于橫波，先出現(xiàn)的波包大于后出現(xiàn)的;而采用均勻化模型只能計(jì)算出一個(gè)波包。采用均勻化材料模型與采用分層離散模型計(jì)算得到的波幅值和相位較為一致，而與采用連續(xù)材料模型計(jì)算的波形結(jié)果差別較大。

采用不同計(jì)算模型而導(dǎo)致波場(chǎng)、時(shí)域響應(yīng)產(chǎn)生差異的主要原因在于，功能梯度結(jié)構(gòu)的宏觀材料屬性沿板厚度方向線性變化，結(jié)構(gòu)中相當(dāng)于存在很多材料界面，使得結(jié)構(gòu)中同時(shí)出現(xiàn)對(duì)稱和反對(duì)稱模態(tài)的波。均勻化模型假設(shè)材料屬性不變，因此只能得到單一模態(tài)的波;分層離散模型，采用分層均勻化的方法建模材料屬性的空間變化特性，模擬精度取決于結(jié)構(gòu)離散程度;連續(xù)模型在單元內(nèi)、外部同時(shí)考慮材料屬性的連續(xù)變化，能夠準(zhǔn)確復(fù)現(xiàn)功能梯度材料結(jié)構(gòu)中的波傳播行為。

圖6 不同材料模型下t=200μs時(shí)結(jié)構(gòu)的位移場(chǎng)，單位μm。Fig.6 Displacements(unit:μm)field with the three models of material properties at 200μs

圖7 不同模型計(jì)算下監(jiān)測(cè)點(diǎn)x，y方向的位移響應(yīng)Fig.7 Displacements of the receiver under the different models in x and y direction

2.2.3 相速度的頻散

高頻導(dǎo)波在板中傳播會(huì)發(fā)生頻散現(xiàn)象，即波速、頻率會(huì)隨著傳播距離發(fā)生變化。進(jìn)一步地，利用板上表面各個(gè)響應(yīng)提取點(diǎn)(均布的41個(gè)點(diǎn))的前400μs x、y方向的時(shí)域響應(yīng)，去除邊界反射波的影響，采用Zero－Crossing方法［15－16］，可以計(jì)算出相速度在傳播過(guò)程中隨頻率的變化，研究不同材料模型對(duì)導(dǎo)波相速度頻散的影響。

圖8(a)、(b)分別給出了采用三種不同模型求解得到的x方向位移響應(yīng)S0模式下、y方向位移響應(yīng)A0模式下的相速度?？梢钥吹剑恳环N模型求解得到的相速度在一定頻域內(nèi)，大小基本保持不變，沒有出現(xiàn)明顯的頻散現(xiàn)象。模型Ⅰ、Ⅱ相速度比較接近;而模型Ⅲ與模型Ⅰ、Ⅱ相比，差異較大。

從以上相速度結(jié)果可以發(fā)現(xiàn)，首先出現(xiàn)的A0、S0模式下波的相速度頻散現(xiàn)象不明顯;不同模型間的相速度數(shù)值差異，說(shuō)明采用均勻化假設(shè)和有限程度的分層離散求解功能梯度材料結(jié)構(gòu)波傳播的相速度的誤差較大。

圖9中給出了采用模型Ⅱ、Ⅲ求解得到的y方向S0模式波的相速度?？梢钥吹?，一定頻率范圍內(nèi)，每一種模型求解得到的相速度大小隨頻率增大而增大，出現(xiàn)明顯的頻散現(xiàn)象。采用分層離散模型計(jì)算得到的相速度高于采用連續(xù)材料模型計(jì)算得到的值。

圖8 三種模型求解得到x方向S0模式、y方向A0模式的相速度Fig.8 Phase velocity of mode S0 in x direction and mode A0 in y direction solved by the three models

圖9 模型Ⅱ、Ⅲ求解得到y(tǒng)方向A0模式的相速度Fig.9 Phase velocity of mode A0 solved by the modelⅡ、Ⅲ in the y direction

3 結(jié)論

本文推導(dǎo)了一種任意四邊形形狀的二維Gauss－Lobatto－Legendre時(shí)域譜單元，采用連續(xù)材料模型在單元內(nèi)部考慮材料屬性的連續(xù)變化，并將其用于超聲導(dǎo)波在功能梯度板結(jié)構(gòu)中的波傳播分析。分別比較了均勻化模型、分層離散、連續(xù)材料模型對(duì)結(jié)構(gòu)波場(chǎng)、時(shí)域響應(yīng)、相速度頻散行為的影響。主要結(jié)論如下:

(1)與理論解比較，驗(yàn)證了推導(dǎo)單元的有效性。

(2)三種模型比較，均勻化模型無(wú)法準(zhǔn)確描述功能梯度材料結(jié)構(gòu)中的波場(chǎng)行為;采用分層離散模型計(jì)算的波響應(yīng)幅值、相位和相速度均與采用連續(xù)材料模型計(jì)算的結(jié)果有差異。采用連續(xù)材料模型能更好模擬功能材料宏觀材料性質(zhì)空間連續(xù)變化的特征。

(3)功能梯度材料中對(duì)稱模式縱波、反對(duì)稱模式橫波的相速度頻散現(xiàn)象不明顯，對(duì)稱模式橫波的相速度頻散明顯。

［1］仲政，吳林志，陳偉球.功能梯度材料與結(jié)構(gòu)的若干力學(xué)問(wèn)題研究進(jìn)展［J］.力學(xué)進(jìn)展，2010，40(5):528－541.ZHONG Zheng，WU Lin－zhi，CHENG Wei－qiu.Progress in the study on mechanics problems of functionally graded materials and structures［J］.Adv Mech，2010，40(5):528－541.

［2］Sun D，Luo S N.Wave propagation and transient response of a FGM plate under a point impact load based on higher－order shear deformation theory［J］.Composite Structures，2011，93(5):1474－1484.

［3］孫丹，羅松南.四邊固支功能梯度板中波的傳播［J］.振動(dòng)與沖擊，2011，30(4):244－247.SUN Dan， LUO Song－nam. Wave propagation in a rectangular functionally graded material plate with clamped supports［J］.Journal of Vibration and Shock，2011，30(4):244－247.

［4］Lefebvre J E，Zhang V，Gazalet J，et al.Acoustic wave propagation in continuous functionally graded plates:an extension of the Legendre polynomial approach［J］.Ultrasonics， Ferroelectrics and Frequency Control， IEEE Transactions on，2001，48(5):1332－1340.

［5］Cao X，Jin F，Jeon I.Calculation of propagation properties of Lamb waves in a functionally graded material(FGM)plate by power series technique［J］.NDT ＆ E International，2011，44(1):84－92.

［6］Patera A T.A spectral element method for fluid dynamics:laminar flow in a channel expansion［J］. Journal of Computational Physics，1984，54(3):468 －488.

［7］Palacz M，Krawczuk M.Analysis of longitudinal wave propagation in a cracked rod by the spectral element method［J］.Computers＆ Structures，2002，80(24):1809－1816.

［8］Kudela P，Krawczuk M，Ostachowicz W.Wave propagation modelling in 1D Structures using spectral finite elements［J］.Journal of Sound and Vibration，2007，300(1):88－100.

［9］Komatitsch D，Martin R，Tromp J，et al.Wave propagation in 2－D elastic media using a spectral element method with triangles and quadrangles［J］.Journal of Computational Acoustics，2001，9(2):703－718.

［10］Peng H，Meng G，Li F.Modeling of wave propagation in plate structures using three－dimensional spectral element method for damage detection［J］.Journal of Sound and Vibration，2009，320(4):942－954.

［11］Chakraborty A，Gopalakrishnan S.A spectrally formulated finite element for wave propagation analysis in functionally graded beams［J］.International Journal of Solids and Structures，2003，40(10):2421 －2448.

［12］Chakraborty A，Gopalakrishnan S.A higher－order spectral element for wave propagation analysis in functionallygraded materials［J］.Acta Mechanica，2004，172(1 －2):17 －43.

［13］Hedayatrasa S，Bui T Q， Zhang C， et al. Numerical modeling of wave propagation in functionally graded materials using time－domain spectral Chebyshev elements［J］.Journal of Computational Physics，2014，258:381－404.

［14］Chen W Q，Wang H M，Bao R H.On calculating dispersion curves of waves in a functionally graded elastic plate［J］.Composite Structures，2007，81(2):233 －242.

［15］Ma?eika L，Draudvilien L，?ukauskas E.Influence of the dispersion on measurement of phase and group velocities of Lamb waves［J］.Ultrasound，2009，64(4):18 －21.

［16］Ma?eika L，Draudvilien L.Analysis of the zero－crossing technique in relation to measurements of phase velocities of the Lamb waves［J］.Ultrasound，2010，66(2):7 －12.