板中熱彈波傳播:一種改進的勒讓德多項式方法1)

王現輝 李方琳 劉宇建 陳會濤 禹建功

(河南理工大學機械與動力工程學院,河南焦作 454000)

引言

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

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

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

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

1 改進的勒讓德多項式法

1.1 分數階熱彈板導波傳播問題求解

基于分數階LS 熱彈理論,熱彈導波傳播問題的控制方程可以如下表達

應變和位移關系,以及本構方程表達式如下

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

為求解該問題,必須進行無量綱化[30]

進一步,假設板內為x+方向傳播的自由諧波,則位移和溫度的表達式可設為

其中,k為波數,ω 是角頻率.在此基礎上,將式(2)～式(5)代入式(1),簡化推導過程,并重寫為,控制方程(1)由可轉化為下面的方程組

在式(6) 中,·′和·′′表示對z的一階和二階導數.為求解方程組(6),可將位移、溫度等物理場變量進行勒讓德多項式近似

式中,Pm,分別是勒讓德多項式及其展開系數,h為板的厚度.在等溫邊界條件中,溫度表達式可表達如下[11]

不同于等溫邊界條件,絕熱邊界條件無法通過一個展開公式進行表達,需再假設溫度梯度的表達式如下

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

通過溫度及其梯度的表達式,可擴展勒讓德多項式法求解絕熱邊界條件下的熱彈問題.

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

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

然而,方程(9) 并不具備一般特征值結構形式.本文弓入波數相關特征向量,進行矩陣變換得到特征值結構形式的方程組.設向量q=kp,方程(9)可轉化為如下特征方程

通過特征值方程組(10),分數階熱彈導波問題的特征值和特征向量即可快速得到.

1.2 解析積分

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

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

為求解I2和I3,參考文獻[31],推導勒讓德多項式導數的展開式,簡化推導過程,其結果是

基于這兩個表達式,當m＞n,mod(m-n,2)=1,l=(m-n-1)/2,I2=2(2m-4l-1)/(2n+1),其他I2=0.當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 個解析積分表達式,可降低算法的計算量,有效減少計算時間,提高計算效率.

2 數值算例

本節將提供3 個算例,第一個算例驗證算法和程序的有效性,第二個算例表明改進算法的計算效率,第三個算例研究分數階次的影響.本節中,板的材料性質和文獻[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).本文數值實驗所使用軟件為Mathematica,電腦配置CPU:Inter core I7-4790,3.6GHz;內存:16G.

2.1 改進算法的有效性

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

2.2 改進算法的高效性

圖1 和已有結果對比Fig.1 Comparing with the existed result

表1 390 頻率計算點的CPU 計算時間(CLPA)Table 1 CPU time with 390 frequencies(CLPA)

表2 390 頻率計算點的CPU 計算時間(AILPA)Table 2 CPU time with 390 frequencies(AILPA)

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

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

2.3 分數階的影響

在本節中,將通過對不同分數階次的頻散曲線和應力、位移、溫度分布進行對比,研究分數階次的影響,松弛時間τ0=1.440×10-13s (無量綱t0=1).圖3 為不同分數階次時的相速度曲線,從圖中可知,熱波波速大于彈波波速.另外,大部分的相速度曲線(彈波相速度)在不同的分數階次下幾乎完全相同,這表明彈波傳播受分數階次的影響較小,這主要因為在熱彈控制方程中,熱彈耦合系數較小,僅為0.002 49,因此溫度對彈波的影響也非常小.同時,圖3 標注的熱波曲線隨著分數階有一定的變化,這表明熱波傳播受分數階的影響較大.但是兩個熱波相速度在分數階的影響下,變化趨勢并不相同.對于無截止頻率的熱波,分數階α 越小,其相速度越大.而對于另一熱波,當頻率較小時,分數階α 越大,相速度越大;當頻率較大時,分數階α 越大,相速度越小.

圖3 不同分數階次的相速度曲線Fig.3 Phase velocity comparison of different fractional order

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

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

圖4 不同分數階的彈波虛部波數對比Fig.4 Comparison of imaginary part wave number of elastic wave with different fractional order

3 結論

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

(1)分數階α 對熱波相速度具有較大的影響,但對彈波的相速度影響較小.然而,分數階對彈波衰減具有一定的影響;

(2)分數階α 影響所有模態的溫度分布,但僅熱波的位移和應力受分數階α 的影響較大.

圖6 應力分布Ω=1Fig.6 Stress distribution with Ω=1

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