王現(xiàn)輝 李方琳 劉宇建 陳會(huì)濤 禹建功

(河南理工大學(xué)機(jī)械與動(dòng)力工程學(xué)院,河南焦作 454000)

引言

超聲導(dǎo)波檢測(cè)技術(shù)具有高效、低成本的特點(diǎn),成為無(wú)損檢測(cè)領(lǐng)域快速發(fā)展的方向之一[1-3].作為一種有效的導(dǎo)波傳播求解方法,勒讓德正交多項(xiàng)式方法直接將邊界條件通過(guò)矩形窗函數(shù)弓入控制方程,簡(jiǎn)化了邊界條件的施加方式; 將求解問(wèn)題的微分方程轉(zhuǎn)化為特征值方程,直接求得表征導(dǎo)波傳播和衰減的復(fù)特征值,具有精確、高效、適用于處理復(fù)雜結(jié)構(gòu)問(wèn)題等優(yōu)點(diǎn).該方法自1999 年提出之后[4],快速應(yīng)用于各種材料結(jié)構(gòu),如多層材料[4-5]、圓柱結(jié)構(gòu)[6-8]、梯度材料[9-10]、熱彈梯度材料[11]、磁電彈耦合問(wèn)題[12]、壓電材料[13]等.在勒讓德正交多項(xiàng)式方法的求解過(guò)程中,需要使用勒讓德多項(xiàng)式近似位移、溫度等物理場(chǎng)變量.同時(shí),為使該方法可解,需要進(jìn)一步讓控制方程推導(dǎo)得到的等式左乘相應(yīng)的勒讓德多項(xiàng)式函數(shù),并對(duì)其從結(jié)構(gòu)的一邊界向另一邊界進(jìn)行積分.由于該積分核函數(shù)中包含有勒讓德多項(xiàng)式及其導(dǎo)數(shù),使得該積分的計(jì)算時(shí)間和計(jì)算量均較大.為克服該缺陷,本文基于勒讓德多項(xiàng)式的正交性,推導(dǎo)所有積分的解析表達(dá)式,降低計(jì)算量,減少計(jì)算時(shí)間,提高計(jì)算效率.

隨著導(dǎo)波無(wú)損檢測(cè)技術(shù)向在線檢測(cè)方向發(fā)展,高溫環(huán)境下的導(dǎo)波技術(shù)受到了越來(lái)越多的重視,高溫下的熱彈耦合問(wèn)題也成為必須要考慮的問(wèn)題.一些廣義熱彈理論,如LS 理論[14],GL 理論[15],GN 理論[16]等被用來(lái)處理此類問(wèn)題.當(dāng)采用勒讓德正交多項(xiàng)式方法求解這類問(wèn)題時(shí),溫度變量需要采用某種形式的勒讓德多項(xiàng)式近似表達(dá).對(duì)于等熱邊界條件,僅需在常規(guī)勒讓德多項(xiàng)式近似的基礎(chǔ)上,乘以z(z-h)項(xiàng)即可[11].然而,對(duì)于其他邊界條件的熱彈問(wèn)題,如絕熱邊界條件問(wèn)題,由于缺乏合適的邊界條件處理方法,勒讓德正交多項(xiàng)式方法無(wú)能為力.為進(jìn)一步擴(kuò)展該方法的應(yīng)用范圍,本文借助于矩形窗函數(shù),發(fā)展一種絕熱邊界條件處理方法.

LS 理論作為一種有效的廣義熱彈理論,近年來(lái)得到了廣泛的發(fā)展[17-19].然而,對(duì)于很多材料和物理過(guò)程,如熱黏彈材料和低溫過(guò)程等,經(jīng)典的整數(shù)階理論難以有效進(jìn)行求解[14].幸運(yùn)的是,分?jǐn)?shù)階熱彈理論[20-23]的發(fā)展為這些問(wèn)題的解決提供了方案.在分?jǐn)?shù)階理論的基礎(chǔ)上,Tiwari 等[24]、Deswal等[25]、Kumar[26-27]對(duì)于平面波傳播問(wèn)題進(jìn)行了研究;許光映等[28]研究了短脈沖激光加熱的溫度場(chǎng)及熱應(yīng)力場(chǎng)的熱物理行為;何天虎等[29]基于非局部效應(yīng)和記記依賴微分修正的廣義熱彈性理論,研究了兩端固定、受移動(dòng)熱源作用的有限長(zhǎng)熱彈桿的動(dòng)態(tài)響應(yīng).然而就作者所知,分?jǐn)?shù)階熱彈理論目前并未應(yīng)用于熱彈導(dǎo)波傳播.在本文中,作者初步將分?jǐn)?shù)階熱彈理論弓入勒讓德正交多項(xiàng)式方法,為以后復(fù)雜材料或者物理過(guò)程的熱彈問(wèn)題有效求解打下基礎(chǔ).

綜上所述,本文提出一種改進(jìn)的勒讓德多項(xiàng)式方法,求解分?jǐn)?shù)階熱彈板中的導(dǎo)波傳播.推導(dǎo)求解方法中積分的解析表達(dá)式,提高計(jì)算效率;弓入溫度梯度表達(dá)式,發(fā)展適合勒讓德多項(xiàng)式級(jí)數(shù)的絕熱邊界條件處理方法.與已有文獻(xiàn)結(jié)果對(duì)比表明改進(jìn)方法的正確性;與已有方法的CPU 計(jì)算時(shí)間對(duì)比說(shuō)明改進(jìn)方法的高效性.最后將改進(jìn)的方法用于求解分?jǐn)?shù)階熱彈板中的導(dǎo)波傳播,研究分?jǐn)?shù)階次對(duì)頻散、衰減曲線和應(yīng)力、位移、溫度分布等的影響.

1 改進(jìn)的勒讓德多項(xiàng)式法

1.1 分?jǐn)?shù)階熱彈板導(dǎo)波傳播問(wèn)題求解

基于分?jǐn)?shù)階LS 熱彈理論,熱彈導(dǎo)波傳播問(wèn)題的控制方程可以如下表達(dá)

應(yīng)變和位移關(guān)系,以及本構(gòu)方程表達(dá)式如下

其中,u是位移,T是溫度,Tij為應(yīng)力,εij是應(yīng)變,Kj是材料常數(shù),βi為熱脹系數(shù),T0=296 K,Ce為特定常數(shù),ρ 是密度,Cij是彈性系數(shù),α 是分?jǐn)?shù)階次,τ0是松弛時(shí)間.

為求解該問(wèn)題,必須進(jìn)行無(wú)量綱化[30]

進(jìn)一步,假設(shè)板內(nèi)為x+方向傳播的自由諧波,則位移和溫度的表達(dá)式可設(shè)為

其中,k為波數(shù),ω 是角頻率.在此基礎(chǔ)上,將式(2)～式(5)代入式(1),簡(jiǎn)化推導(dǎo)過(guò)程,并重寫(xiě)為,控制方程(1)由可轉(zhuǎn)化為下面的方程組

在式(6) 中,·′和·′′表示對(duì)z的一階和二階導(dǎo)數(shù).為求解方程組(6),可將位移、溫度等物理場(chǎng)變量進(jìn)行勒讓德多項(xiàng)式近似

式中,Pm,分別是勒讓德多項(xiàng)式及其展開(kāi)系數(shù),h為板的厚度.在等溫邊界條件中,溫度表達(dá)式可表達(dá)如下[11]

不同于等溫邊界條件,絕熱邊界條件無(wú)法通過(guò)一個(gè)展開(kāi)公式進(jìn)行表達(dá),需再假設(shè)溫度梯度的表達(dá)式如下

將式(7)中的溫度表達(dá)式X(z)對(duì)z向求導(dǎo),其結(jié)果必然要等于式(8)中的溫度梯度表達(dá)式,然后讓該等式兩邊同乘以Qj(z),j=0,1,...,N,并進(jìn)行積分,則可以得到方程組Hp5=Hp4.易知,p5=p4.

通過(guò)溫度及其梯度的表達(dá)式,可擴(kuò)展勒讓德多項(xiàng)式法求解絕熱邊界條件下的熱彈問(wèn)題.

將式(7)和式(8)弓入方程(6),方程可轉(zhuǎn)為下面的矩陣形式

Ai j,Bi j,Cij,Mij,D,E是將方程(7)和(8)代入式(6),并進(jìn)行積分后形成的矩陣.方程(9a)由方程(6a),(6c) 和(6d) 得到,方程(9b) 由獨(dú)立方程(6b) 推導(dǎo)求得.

然而,方程(9) 并不具備一般特征值結(jié)構(gòu)形式.本文弓入波數(shù)相關(guān)特征向量,進(jìn)行矩陣變換得到特征值結(jié)構(gòu)形式的方程組.設(shè)向量q=kp,方程(9)可轉(zhuǎn)化為如下特征方程

通過(guò)特征值方程組(10),分?jǐn)?shù)階熱彈導(dǎo)波問(wèn)題的特征值和特征向量即可快速得到.

1.2 解析積分

在程序運(yùn)行過(guò)程中,積分所消耗的時(shí)間占總CPU 計(jì)算時(shí)間的比例可高達(dá)97%以上(見(jiàn)表2),因此有必要通過(guò)有效的方法,降低積分計(jì)算量,減少積分時(shí)間,提高計(jì)算效率.積分解析化是實(shí)現(xiàn)該目標(biāo)的一種有效手段.本節(jié)基于勒讓德多項(xiàng)式的正交性,對(duì)程序運(yùn)行過(guò)程中的所有積分進(jìn)行解析表達(dá).通過(guò)分析,所有積分可歸納為以下5 種積分

h(t)為Heaviside 函數(shù).在這5 種積分中,I1,I4和I5可直接得到

為求解I2和I3,參考文獻(xiàn)[31],推導(dǎo)勒讓德多項(xiàng)式導(dǎo)數(shù)的展開(kāi)式,簡(jiǎn)化推導(dǎo)過(guò)程,其結(jié)果是

基于這兩個(gè)表達(dá)式,當(dāng)m＞n,mod(m-n,2)=1,l=(m-n-1)/2,I2=2(2m-4l-1)/(2n+1),其他I2=0.當(dāng)m＞n+1,mod(m-n,2)=0,p=(m-n-2)/2,I3=2(m-4p-3)(p+1)(2m-2p-1)/(2n+1),其他I3=0.

基于這5 個(gè)解析積分表達(dá)式,可降低算法的計(jì)算量,有效減少計(jì)算時(shí)間,提高計(jì)算效率.

2 數(shù)值算例

本節(jié)將提供3 個(gè)算例,第一個(gè)算例驗(yàn)證算法和程序的有效性,第二個(gè)算例表明改進(jìn)算法的計(jì)算效率,第三個(gè)算例研究分?jǐn)?shù)階次的影響.本節(jié)中,板的材料性質(zhì)和文獻(xiàn)[30] 一致,C11,C13,C33,C44,C55分別為5.74×1011,1.27×1011,4.33×1011,1.19×1011,1.08×1011N/m2,P為3.2×103kg/m3,Ce為670 J·kg·(°)/m,β1和β3分別為3.22×106和2.71×106N/((°)·m2),k1和k2分別為55.4 和43.5 W/(m·K).本文數(shù)值實(shí)驗(yàn)所使用軟件為Mathematica,電腦配置CPU:Inter core I7-4790,3.6GHz;內(nèi)存:16G.

2.1 改進(jìn)算法的有效性

由于缺乏分?jǐn)?shù)階熱彈板中導(dǎo)波傳播的研究成果,本文將和整數(shù)階(α=1)的熱彈導(dǎo)波頻散曲線[30]進(jìn)行了對(duì)比,如圖1 所示.結(jié)果顯示由改進(jìn)的勒讓德多項(xiàng)式法得到的頻散曲線和已有結(jié)果完全一致,表明本文的算法和程序是正確有效的.

2.2 改進(jìn)算法的高效性

圖1 和已有結(jié)果對(duì)比Fig.1 Comparing with the existed result

表1 390 頻率計(jì)算點(diǎn)的CPU 計(jì)算時(shí)間(CLPA)Table 1 CPU time with 390 frequencies(CLPA)

表2 390 頻率計(jì)算點(diǎn)的CPU 計(jì)算時(shí)間(AILPA)Table 2 CPU time with 390 frequencies(AILPA)

在本節(jié)中,通過(guò)對(duì)不同展開(kāi)階次勒讓德多項(xiàng)式的CPU 計(jì)算時(shí)間進(jìn)行統(tǒng)計(jì),研究改進(jìn)方法的計(jì)算效率.傳統(tǒng)勒讓德正交多項(xiàng)式方法(conventional Legendre polynomial approach,CLPA)和本文改進(jìn)的勒讓德正交多項(xiàng)式方法(analytical integration Legendre polynomial approach,AILPA) CPU 計(jì)算時(shí)間(s)分別如表1 和表2 所示.表1 結(jié)果表明,已有方法積分時(shí)間占總計(jì)算時(shí)間的比重最高為97.8%,這意味著幾乎程序所有的計(jì)算時(shí)間都用來(lái)進(jìn)行積分.表2 結(jié)果表明積分時(shí)間占總CPU 計(jì)算時(shí)間的比重最高為17.2%,這意味著僅有不足1/5 的計(jì)算時(shí)間都用來(lái)進(jìn)行積分.表1 和表2 的數(shù)據(jù)對(duì)比表明改進(jìn)方法計(jì)算積分的時(shí)間大大減少,致使總體CPU 計(jì)算時(shí)間大幅度降低.圖2 為AILPA 的總計(jì)算時(shí)間占CLPA 的總計(jì)算時(shí)間的百分比,從圖上可知,當(dāng)N=5 時(shí),占比最高,此時(shí)為21.05%,并且隨著展開(kāi)階N的增大,該占比不斷下降.值得注意的是,在N=20 時(shí),AILPA 的總計(jì)算時(shí)間僅為10.2 s,為CLPA 總計(jì)算時(shí)間(392.4 s) 的2.6%.因此,和CLPA 相比,改進(jìn)算法的計(jì)算效率得到極大提高.

圖2 兩種方法CPU 總時(shí)間的比值Fig.2 Ratio of two total CPU time

2.3 分?jǐn)?shù)階的影響

在本節(jié)中,將通過(guò)對(duì)不同分?jǐn)?shù)階次的頻散曲線和應(yīng)力、位移、溫度分布進(jìn)行對(duì)比,研究分?jǐn)?shù)階次的影響,松弛時(shí)間τ0=1.440×10-13s (無(wú)量綱t0=1).圖3 為不同分?jǐn)?shù)階次時(shí)的相速度曲線,從圖中可知,熱波波速大于彈波波速.另外,大部分的相速度曲線(彈波相速度)在不同的分?jǐn)?shù)階次下幾乎完全相同,這表明彈波傳播受分?jǐn)?shù)階次的影響較小,這主要因?yàn)樵跓釓椏刂品匠讨?熱彈耦合系數(shù)較小,僅為0.002 49,因此溫度對(duì)彈波的影響也非常小.同時(shí),圖3 標(biāo)注的熱波曲線隨著分?jǐn)?shù)階有一定的變化,這表明熱波傳播受分?jǐn)?shù)階的影響較大.但是兩個(gè)熱波相速度在分?jǐn)?shù)階的影響下,變化趨勢(shì)并不相同.對(duì)于無(wú)截止頻率的熱波,分?jǐn)?shù)階α 越小,其相速度越大.而對(duì)于另一熱波,當(dāng)頻率較小時(shí),分?jǐn)?shù)階α 越大,相速度越大;當(dāng)頻率較大時(shí),分?jǐn)?shù)階α 越大,相速度越小.

圖3 不同分?jǐn)?shù)階次的相速度曲線Fig.3 Phase velocity comparison of different fractional order

由于彈性模態(tài)更適合于進(jìn)行無(wú)損檢測(cè),因此需進(jìn)一步研究分?jǐn)?shù)階次對(duì)彈波衰減的影響.圖4 為不同分?jǐn)?shù)階次時(shí)彈波波數(shù)的虛部值對(duì)比.結(jié)果顯示,隨著分?jǐn)?shù)階次的變化,彈波波數(shù)的虛部值也隨著變化.其中模態(tài)2 和3 受分?jǐn)?shù)階的影響較大,且具有相同的趨勢(shì),即分?jǐn)?shù)階α 越大,波數(shù)虛部值越大.這意味較大的分?jǐn)?shù)階α 值,在部分彈波模態(tài)傳播時(shí)具有較大的衰減速度.因此,在使用彈波進(jìn)行無(wú)損檢測(cè)時(shí),分?jǐn)?shù)階對(duì)衰減的影響不可忽略.

圖5～圖7 分別給出了位移(uj)、應(yīng)力(Tij) 和溫度(T)及溫度梯度(Tz)分布的情況,其中—–表示α=0.2,–·–·–表示α=0.5,--------表示α=0.8.從3幅圖上可知,(1)分?jǐn)?shù)階α 影響位移、應(yīng)力和溫度的分布; (2) 分?jǐn)?shù)階α 影響所有模態(tài)的溫度分布,但僅有熱波的位移和應(yīng)力受到分?jǐn)?shù)階α 的影響較大;(3)z向應(yīng)力和溫度梯度分布曲線均是從0 到0,這進(jìn)一步表明了程序的正確性.另外,溫度梯度分布曲線從0到0 說(shuō)明絕熱邊界條件處理方法的正確性.

圖4 不同分?jǐn)?shù)階的彈波虛部波數(shù)對(duì)比Fig.4 Comparison of imaginary part wave number of elastic wave with different fractional order

3 結(jié)論

針對(duì)勒讓德正交多項(xiàng)式方法求解熱彈導(dǎo)波傳播時(shí)存在的兩個(gè)不足之處,本文提出一種改進(jìn)的勒讓德正交多項(xiàng)式方法,以求解分?jǐn)?shù)階熱彈板中的導(dǎo)波傳播.本文主要?jiǎng)?chuàng)新點(diǎn)是:(1) 采用解析積分代替原有的數(shù)值積分,極大提高了算法的計(jì)算效率.該積分表達(dá)式可推廣至壓電材料、圓柱等均勻結(jié)構(gòu); (2)發(fā)展一種適用于勒讓德多項(xiàng)式級(jí)數(shù)的絕熱邊界條件處理方法,擴(kuò)展了原有算法的求解范圍;(3)將改進(jìn)的方法應(yīng)用于分?jǐn)?shù)階熱彈板,研究分?jǐn)?shù)階次對(duì)頻散、衰減曲線和位移、應(yīng)力、溫度分布等的影響.主要研究結(jié)論如下.

(1)分?jǐn)?shù)階α 對(duì)熱波相速度具有較大的影響,但對(duì)彈波的相速度影響較小.然而,分?jǐn)?shù)階對(duì)彈波衰減具有一定的影響;

(2)分?jǐn)?shù)階α 影響所有模態(tài)的溫度分布,但僅熱波的位移和應(yīng)力受分?jǐn)?shù)階α 的影響較大.

圖6 應(yīng)力分布Ω=1Fig.6 Stress distribution with Ω=1

圖7 溫度及其梯度分布Ω=1Fig.7 Temperature and its gradient distribution with Ω=1