基于時域譜單元的功能梯度材料結構波傳播分析

徐 超，王 騰

(西北工業大學 航天學院，西安 710072)

功能梯度材料的宏觀材料特性表現出空間連續變化的性質。材料成分和性質的變化會改變彈性波在結構中的傳播行為，研究彈性波在功能梯度等非均勻材料結構中的傳播特性對發展新材料結構的損傷檢測和健康監測方法具有重要的意義［1］。

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

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

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

1 任意四邊形譜單元推導

1.1 位移、形狀插值

如圖1所示，與傳統有限元不同，自然坐標系下標準譜單元的節點非等距地分布于區域Λ∈［－1，1］2內，單元上節點坐標位置根據特殊多項式的零點確定。非等距節點布置能有效克服等距節點布置在插值時引起的龍格效應。這里采用 Gauss－Lobatto－Legendre多項式確定單元中各節點的坐標位置。以ξ方向為例，Gauss－Lobatto－Legendre 多項式為

式中:Pn－1(ξ)為n－1階勒讓德多項式。單元階次不同，對應的節點坐標位置不同，如當n=5時，節點ξ方向坐標為

圖1 2D譜單元及5×5節點分布Fig.1 2D spectral element with 5 ×5 nodes

同理，單元η方向上的節點坐標也可類似確定。

仍采用拉格朗日插值函數對單元內位移場進行插值，即

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

實際應用中，剖分離散后的單元幾何形狀可能不是標準的邊長為1的正方形。考慮單元形狀為任意四邊形的情況，需進行形狀插值建立任意四邊形單元與標準譜單元的幾何映射關系。

設任意四邊形四個角節點的坐標為(xk，yk)，k=1，2，3，4，采用一階拉格朗日插值函數建立單元形狀插值關系:

1.2 功能梯度材料模型

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

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

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

(3)連續材料模型:將材料屬性考慮為坐標的連續函數，則單元內材料性質也是坐標的連續函數，材料密度矩陣μe和彈性矩陣De可表示為:

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

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

式中:yu，yd代表結構上下界面的坐標，h=yu－yd。χ表示材料分布的不同梯度，當χ=0為均勻材料。

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

1.3 單元剛度和質量矩陣

利用上述位移和形狀插值函數，根據Hamilton原理建立變分方程，可推導獲得單元質量矩陣Me，單元剛度矩陣Ke和單元上的等效外力Fe表達式為:

式中:Be矩陣為幾何矩陣，可以由位移－應變幾何關系求出:

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

式(7)中J為雅克比矩陣，描述任意四邊形單元與標準譜單元的幾何關系，根據式(3)有

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

由于勒讓德多項式的正交性，單元質量矩陣為對角矩陣，這可顯著降低波動力響應時域求解的計算耗費，較切比雪夫多項式時域譜單元有明顯優勢。

2 數值算例

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

為了驗證所推導單元的有效性，考慮一中心受載的無限大FGM板。文獻［2］采用彈性板殼理論，考慮高階剪切變形與轉動慣量的影響，推導獲得了極坐標下的解析解。本文采用平面應變模型，利用前文推導的連續材料模型譜單元離散板結構，將獲得的數值解與解析解進行對比。

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

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

應用譜單元數值求解時，單元尺寸為0.01×0.067 m，每個單元節點數為5×5，積分步長取2×10－8s，監測距加載處x=2 m處的響應和解析解對比。采用文獻［2］中的無量綱方法，無量綱后的時域響應如圖3(a)、(b)所示。其中，圖3(a)給出了x方向的位移響應，圖3(b)給出了y方向的位移響應。由圖可知，數值解和解析解吻合較好，驗證了本文推導單元的有效性。需要說明的是，這里是依據平面應變假設進行數值求解，給出了近似解，目的是驗證單元的波傳播描述能力。由于連續材料模型的單元精度也受單元規模和積分階次制約，給出的解也是數值近似解，因此圖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 面內一端受水平激勵的FGM板

2.2.1 模型描述

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

圖4 板xy截面及計算網格Fig.4 Plane xy of plate and computational grids

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

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

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

2.2.2 波場及時域響應

左端施加激勵后，導波會沿著板截面傳播，同時受到上下邊界的約束。

圖6給出了采用三種不同的材料模型計算出的結構在200μs時x、y方向的位移場云圖。在二維結構中，由于泊松效應，水平方向的激勵在水平和垂直方向都引起結構的振動。在水平方向，能量主要集中在以對稱模式S0波傳播的模態上;在垂直方向，能量主要集中以反對稱模式A0波傳播的模態上。

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

為了進一步定量研究不同模型對波動響應的影響，取板上表面的中點作為數值結果時域響應的提取點(因一般壓電傳感器貼在結構上表面，監測局部響應)，給出其0～400μs的位移響應，結果如圖7所示。

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

采用不同計算模型而導致波場、時域響應產生差異的主要原因在于，功能梯度結構的宏觀材料屬性沿板厚度方向線性變化，結構中相當于存在很多材料界面，使得結構中同時出現對稱和反對稱模態的波。均勻化模型假設材料屬性不變，因此只能得到單一模態的波;分層離散模型，采用分層均勻化的方法建模材料屬性的空間變化特性，模擬精度取決于結構離散程度;連續模型在單元內、外部同時考慮材料屬性的連續變化，能夠準確復現功能梯度材料結構中的波傳播行為。

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

圖7 不同模型計算下監測點x，y方向的位移響應Fig.7 Displacements of the receiver under the different models in x and y direction

2.2.3 相速度的頻散

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

圖8(a)、(b)分別給出了采用三種不同模型求解得到的x方向位移響應S0模式下、y方向位移響應A0模式下的相速度。可以看到，每一種模型求解得到的相速度在一定頻域內，大小基本保持不變，沒有出現明顯的頻散現象。模型Ⅰ、Ⅱ相速度比較接近;而模型Ⅲ與模型Ⅰ、Ⅱ相比，差異較大。

從以上相速度結果可以發現，首先出現的A0、S0模式下波的相速度頻散現象不明顯;不同模型間的相速度數值差異，說明采用均勻化假設和有限程度的分層離散求解功能梯度材料結構波傳播的相速度的誤差較大。

圖9中給出了采用模型Ⅱ、Ⅲ求解得到的y方向S0模式波的相速度。可以看到，一定頻率范圍內，每一種模型求解得到的相速度大小隨頻率增大而增大，出現明顯的頻散現象。采用分層離散模型計算得到的相速度高于采用連續材料模型計算得到的值。

圖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方向A0模式的相速度Fig.9 Phase velocity of mode A0 solved by the modelⅡ、Ⅲ in the y direction

3 結論

本文推導了一種任意四邊形形狀的二維Gauss－Lobatto－Legendre時域譜單元，采用連續材料模型在單元內部考慮材料屬性的連續變化，并將其用于超聲導波在功能梯度板結構中的波傳播分析。分別比較了均勻化模型、分層離散、連續材料模型對結構波場、時域響應、相速度頻散行為的影響。主要結論如下:

(1)與理論解比較，驗證了推導單元的有效性。

(2)三種模型比較，均勻化模型無法準確描述功能梯度材料結構中的波場行為;采用分層離散模型計算的波響應幅值、相位和相速度均與采用連續材料模型計算的結果有差異。采用連續材料模型能更好模擬功能材料宏觀材料性質空間連續變化的特征。

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

［1］仲政，吳林志，陳偉球.功能梯度材料與結構的若干力學問題研究進展［J］.力學進展，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］.振動與沖擊，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.