陳薈鍵, 朱清鋒, 苗鴻臣, 馮志強
(西南交通大學 力學與航空航天學院, 成都 611756)
超聲導波具有傳播距離遠、能量衰減小、對結構損傷敏感等優(yōu)點,在結構健康監(jiān)測和無損檢測領域被廣泛關注[1].傳統(tǒng)的線性導波技術基于導波的線性效應,主要通過測試聲速、衰減、反射率和透射率等參數(shù)來檢測缺陷.然而,線性導波通常僅能檢測與波長相當?shù)娜毕?如孔洞、宏觀裂紋等,而對材料中的初期損傷不敏感.非線性超聲導波兼具線性導波快速檢測和非線性超聲對材料初始損傷高靈敏度識別的雙重優(yōu)點[2].以接觸聲非線性效應為例,超聲導波與微裂紋作用過程中的接觸聲非線性效應可以引起高次諧波、零頻響應等非線性超聲特征[3],這些特征可以用于識別線性超聲無法識別的微裂紋.
板殼結構中的零階水平剪切波(SH0波)因其獨特的非頻散特性,在結構健康監(jiān)測和無損檢測中具有重要的應用前景[4].洞悉SH0波與裂紋的作用規(guī)律對于發(fā)展基于SH0波的裂紋檢測技術具有重要的意義.當前SH0波與裂紋的線性散射規(guī)律已經得到了較為深入的研究[5-8],而關于SH0波與微裂紋的接觸聲非線性散射規(guī)律的研究成果還比較匱乏,特別是結構外部荷載對SH0波非線性散射特征的影響規(guī)律尚未見報道.因此,有必要研究受載條件下SH0波與微裂紋作用的非線性散射規(guī)律.由于此時結構中存在初始應力場和聲場的耦合,使得SH0波與微裂紋作用時的接觸聲非線性作用十分復雜,因此需要借助數(shù)值仿真對該問題進行分析.有限元法[9]具有易實現(xiàn)、靈活性高等優(yōu)點,是模擬結構中導波傳播最常用的方法之一.然而,由于導波具有頻率高、波長短的特點,使得傳統(tǒng)有限元法在求解此問題時,對網格密度的要求高,計算效率較低.時域譜元法[10]結合了偽譜法高階多項式快速收斂性和有限元法復雜幾何適應性的優(yōu)點,能夠以較小的計算代價求解導波傳播問題.然而,譜單元對復雜幾何結構的離散能力不如有限元法.
為求解接觸聲非線性問題,需要準確描述在超聲波作用下裂紋面的張開、閉合以及滑移等接觸摩擦行為[3].目前求解接觸問題主要的方法有罰函數(shù)法和Lagrange乘子法.罰函數(shù)方法[11-12]容易實現(xiàn)且收斂速度快,但是罰因子的選取較為困難.若選取的罰因子過小,則接觸穿透量較大,從而影響數(shù)值精度;若增大罰因子,雖然可以更加符合接觸條件,但是這樣又會引起局部接觸計算的不穩(wěn)定,同時時間步長還需設置得非常小.此外,罰函數(shù)法難以考慮接觸計算時的多點耦合效應,而在實際的接觸計算中,接觸點之間是相互影響的.Lagrange乘子法[13]能夠精確地滿足物體間的法向接觸邊界,但是系統(tǒng)中由于引入了Lagrange乘子,使得方程變量增多,從而增加了整個問題的計算量.由De Saxcé和Feng[14-15]提出的雙勢接觸方法,在不增加自由度的前提下,既能夠準確地滿足法向的Signorini接觸條件,又能滿足切向的Coulomb摩擦定律,在處理接觸問題上具有較高的精度和效率.Chen等[16]結合譜元法和雙勢接觸方法,發(fā)展了可以高效計算接觸聲非線性問題的雙勢譜方法,證明了雙勢方法可以準確處理波動接觸問題.
本文在雙勢譜方法的基礎上,進一步通過mortar方法將譜單元和有限單元進行耦合,以此研究了SH0波與裂紋作用的非線性散射規(guī)律.首先,介紹了處理有限單元/譜單元耦合界面的mortar方法的基本理論.其次,介紹了用于刻畫導波與微裂紋作用時,裂紋面的接觸摩擦行為的雙勢接觸理論.接著,本文介紹了一種預處理半顯式算法,該算法囊括了包含附加界面力的局部接觸力的隱式求解與全局位移的顯式求解,其可以在保證計算效率的同時而不犧牲數(shù)值精度.最后,通過發(fā)展的數(shù)值方法計算了自由狀態(tài)和單軸受載狀態(tài)下SH0波與微裂紋作用的非線性散射規(guī)律,并分析了應力大小對非線性散射場的影響.
Mortar方法是一種非協(xié)調區(qū)域分解算法,它允許將求解區(qū)域分解為多個子區(qū)域,對各個子域分別進行建模和剖分,然后在各個區(qū)域以最適合子域特征的方式進行離散.在mortar方法中,各個子域交界面處的邊界節(jié)點不需要嚴格逐點匹配,而是通過約束方程弱形式建立具有最佳精度的mortar條件來使得界面處的場變量得以傳遞,以此來保證不同子域之間的連續(xù)性.Bernardi等[17]最先于1994年提出了mortar元法的概念.目前mortar方法已被成功應用于處理彈性動力學問題[18-19]、流體力學問題[20]、接觸摩擦問題[21-22]等.
如圖1所示,考慮由兩個互不重疊的子域ΩFE和ΩSFE構成的域Ω,滿足Ω=ΩFE∪ΩSFE和ΩFE∩ΩSFE=?.其中Ω為研究總域;ΩFE采用有限單元離散,稱為有限元區(qū)域;ΩSFE采用譜單元離散,稱為譜單元區(qū)域.兩個子域之間的交界面用γ表示,即γ=?ΩFE∩?ΩSFE.由于兩個子域采用不同的單元類型進行離散,因此在界面γ處存在不匹配性.這種不匹配性主要體現(xiàn)在兩個方面:一是幾何上的不匹配.由于有限單元和譜單元在網格劃分以及單元類型上的差異,界面處兩個區(qū)域的單元節(jié)點往往是不重合的.二是插值函數(shù)上的不匹配.有限元采用的是低階的插值形函數(shù),而譜單元采用高階多項式作為插值形函數(shù),因此兩者在界面兩側擁有不同的形函數(shù)空間.
圖1 有限元和譜單元耦合示意圖Fig. 1 Schematic diagram of the coupling of finite elements and spectral elements
考慮界面γ處的位移連續(xù)性條件:
(1)
式中,uFE和uSFE分別為有限元區(qū)域和譜單元區(qū)域的位移場;Ni是和界面節(jié)點相關的有限單元的形函數(shù);i表示界面的第i個節(jié)點;Nf為有限元區(qū)域界面節(jié)點的個數(shù).為了使用顯式時間積分,位移連續(xù)性條件用等價的速度連續(xù)性條件替代,即
(2)
上式約束條件一般可稱為mortar條件,關于mortar條件的數(shù)學描述可參考文獻[18,23].
Mortar方法的目標是在整個界面γ上滿足連續(xù)性條件和運動平衡方程.既要保證兩個子域在界面上不會發(fā)生分離和重疊,又要保證整個界面的動力平衡.后者可以通過引入適當?shù)慕缑嫦嗷プ饔昧韺崿F(xiàn),本文采用Lagrange乘子法求解該界面力.在不考慮接觸力和材料阻尼的條件下,系統(tǒng)的控制方程可以寫為[18]
(3)
(4)
其中,G為界面耦合矩陣,其只和界面上節(jié)點的初始位置相關,在整個動力學求解過程只需計算一次[19];Rγ為附加界面力矢量,可通過Lagrange乘子法求解得到
Rγ=GTλ,
(5)
式中,λ為Lagrange乘子,將式(5)代入式(3),并對等式兩邊同時乘以GM-1,可以得到
(6)
引入中心差分格式
(7)
(8)
與中間時刻速度表達式
(9)
上式兩端同時乘以G,并做等式變換可得
(10)
將上式代入式(6),并結合式(4)可得
(11)
至此,Lagrange乘子λ可按照常規(guī)線性方程組Ax=B的格式進行求解.一旦求解得到λ,便可通過式(5)求解得到附加界面力矢量Rγ,然后將求解得到的Rγ代入式(3),再通過時間積分算法便可獲得位移場.
雙勢理論[14]由De Saxcé和Feng于1991年提出,該理論以Legendre定理為基礎,通過Fenchel不等式變換,建立了一組能夠處理對偶變量的方程.通過對材料模型能量勢函數(shù)形式進行分析,將材料分為顯式標準材料和隱式標準材料,在此基礎上,提出了雙勢函數(shù)的概念.
對于接觸問題,雙勢函數(shù)Fb可以寫成如下形式:
(12)
(13)
(14)
基于式(12),可以得到接觸問題中雙勢函數(shù)的次法向準則,即
(15)
νFb(-xα,rα*)-νFb(-xα,rα)≥[(rα-νxα)-rα]·(rα*-rα),
(16)
式中,ν>0,其值的選擇確保了數(shù)值計算的收斂性,一般可由減縮的接觸柔度矩陣的最小特征值決定[25].從上式的右端項可知,rα為增廣Lagrange接觸力rα*的臨近點.因此,可將增廣Lagrange接觸力rα*表示為
(17)
式中,rα*在接觸迭代計算中又被稱為預測接觸力.得到預測接觸力后再在Coulomb摩擦錐上進行投影,從而得到真實接觸力,即
rα=PKμ(rα*).
(18)
圖2 Coulomb摩擦錐示意圖Fig. 2 Schematic diagram of Coulomb’s frictional cone
(19)
該投影方程已被嚴格證明等效于不等式(16),并且在數(shù)值計算中,具有出良好的精確性和魯棒性[26].
在考慮接觸力以及附加界面力的情況下,系統(tǒng)的控制方程可以寫為
(20)
引入有限差分格式:
(21)
將式(21)代入式(20),可以得到
(22)
對于線彈性材料,一般可取C=αM,其中α為質量阻尼系數(shù).代入式(22),可預定義以下參數(shù):
(23)
此時可將式(22)重寫為
ut+Δt=(ψ1+ψ2)-1·[(Fext)t-(Fint)t+(Rγ)t+(Rc)t+Δt+2ψ1ut-(ψ1-ψ2)ut-Δt].
(24)
上式可進一步寫成關于每一個節(jié)點自由度i的顯式表達式:
(25)
從上式可以看出,即便在附加界面力通過Lagrange乘子法求解得到后,全局位移和全局接觸力仍為未知量,不能直接求解.本文采取分開求解的策略,先引入接觸映射關系,通過隱式迭代算法求解局部接觸力r,然后通過局部接觸力與全局接觸力的關系將局部接觸力轉換為全局接觸力Rc,最后將全局接觸力引入式(24)中,即可顯式地求解剩余的未知量u.
在接觸運動學中,局部接觸力r與全局接觸力Rc以及局部接觸位移x與全局接觸位移u之間的關系可寫為[16]
(26)
式中,g為接觸系統(tǒng)中的接觸點與投影點之間的初始間隙矢量.
結合式(24)和(26),有
x=T(ψ1+ψ2)-1TTrt+Δt+
T(ψ1+ψ2)-1[(Fext)t-(Fint)t+(Rγ)t+2ψ1ut-(ψ1-ψ2)ut-Δt]+g.
(27)
定義
(28)
式(27)可以進一步寫為
(29)
其中,W是接觸系統(tǒng)中基于質量矩陣構建的Delassus算子.需要注意的是,由于質量矩陣的對角性質,W的計算將會非常迅速.在包含Nc個接觸點的減縮接觸系統(tǒng)中,對每一個接觸點α而言,局部相對位移xα和局部接觸力rα之間有如下關系:
(30)
式中,Wαβ=Tα(ψ1+ψ2)-1(Tβ)T為影響矩陣,其描述了接觸點α和β之間的接觸耦合效應.換言之,對于每一個接觸點α,都考慮了其他接觸點β(β≠α)對α的接觸力的貢獻.
結合式(18)和(30),減縮系統(tǒng)下接觸力的求解可以轉換為尋找χ以滿足f(χ)=0,即
(31)
其中
(32)
(33)
隱式方程(32)可采用Uzawa算法[27]進行局部迭代求解,該算法主要包括“預測-修正”步驟,即
(34)
式中,i和i+1為迭代數(shù).如前所述,通過Uzawa算法求得局部接觸力,然后通過局部接觸力與全局接觸力的關系式(26)將局部接觸力轉換為全局接觸力,最后將全局接觸力代入式(25)中,即可顯式地求解未知量u.
圖3為SH0波與不同角度微裂紋作用的非線性散射場的計算模型.板的面內尺寸為100 mm×100 mm,采用平面應力模型進行計算,所用材料為鋁,其彈性模量為70 GPa,Poisson比為0.3,質量密度為2 700 kg/m3.在板的上下邊界施加預拉/壓應力p,同時在板的四周布置長度為10 mm的漸增阻尼邊界吸收層[28],最大阻尼系數(shù)αmax=107.坐標原點設置在板中心位置.微裂紋角度θ為裂紋長度方向與x軸正方向的夾角,微裂紋幾何中心與板中心重合,其長度設置為6 mm.采用有限單元/譜單元混合單元對板進行離散,其中含裂紋的圓形區(qū)域(半徑R=20 mm)采用有限單元離散,不含缺陷的外部區(qū)域采用譜單元離散.入射SH0波由板左邊界長60 mm的對稱段上施加沿y方向的激勵力產生,激勵波形為Hanning窗調制的五周期正弦波,中心頻率為500 kHz,幅值為100 N.為方便提取散射場信號,取有限元區(qū)域的外邊界節(jié)點作為散射場的響應點.
圖3 SH0波與微裂紋作用的非線性散射場計算模型Fig. 3 The calculation model for the nonlinear scattering characteristics of the SH0 wave encountering cracks in prestressed plates
如圖4所示,得益于本文介紹的處理有限元/譜單元耦合界面的mortar方法,當研究裂紋角度對SH0波散射特性的影響時,可以把外部譜單元區(qū)域固定,然后將內部包含微裂紋的有限元區(qū)域進行整體旋轉,即可得到新的裂紋模型,旋轉后該區(qū)域的節(jié)點坐標可通過初始節(jié)點坐標向量做簡單的張量變換求得.該策略只需要建立一次裂紋模型便可實現(xiàn)任意角度的裂紋建模,因而可以巧妙地避免傳統(tǒng)方法在處理此類問題時對裂紋進行重復建模、網格重新劃分的麻煩.
圖4 不同角度裂紋建模原理圖Fig. 4 Schematic diagram of the crack modeling at different angles
本小節(jié)探究了自由板(p=0)中SH0波與微裂紋作用的非線性散射規(guī)律.取θ=60°,圖5展示了SH0波與微裂紋作用前后的總位移云圖.可以觀察到:① 有限元/譜單元交界面處網格不匹配對SH0波在結構中的傳播沒有任何影響,驗證了不同單元類型區(qū)域之間的連續(xù)性,證明了本文中mortar方法的有效性;② 當SH0波傳到邊界時,絕大部分能量被吸收掉,反射信號基本消失,證明了漸增阻尼邊界吸收層的有效性;③ 當SH0波與微裂紋作用后,會向外輻射一個新的散射聲場,該聲場在不同方向有不同的信號強度.
圖5 不同時刻的非線性散射總位移云圖Fig. 5 The nonlinear scattering fields at different moments
圖6為提取含微裂紋結構中M點的時域信號以及對應的頻譜圖.作為對比,圖6中同時給出了對應的不含裂紋即完整結構的計算結果.從圖6(a)可以看出,兩組信號僅在幅值上有所不同,但波形上的差異并不明顯.為進一步區(qū)分這兩種情況,我們將兩組時域信號做快速Fourier變換,可以得到對應的頻譜圖,如圖6(b)所示.從頻譜圖可以直觀地分辨出兩種工況存在不同的頻率分量.具體地,對于無損傷結構,頻譜圖中只包含500 kHz激勵頻率分量;而對于含微裂紋的結構,在頻譜圖中,除激勵頻率500 kHz之外,還產生了1 000 kHz,1 500 kHz,2 000 kHz等明顯的高次諧波分量.該結果表明,可以通過SH0波與微裂紋相互作用產生的高次諧波分量來對微裂紋進行有效識別.
圖6 無損傷和含微裂紋損傷工況下,時域信號和對應頻譜圖的比較Fig. 6 Comparison of time domain signals and corresponding frequency spectrums between the pristine case and the microcrack damaged case
實際情況下的微裂紋的取向往往是未知的,而工程中對裂紋擴展方向的判斷又至關重要[29].因此,建立非線性散射規(guī)律與微裂紋取向角之間的關系,對于裂紋檢測傳感器的布置以及傳感信號識別具有重要的意義.需要注意的是,由于SH0波僅有一個垂直于傳播方向的面內位移分量,因此極坐標系下的切向位移可以用來表示SH0波模態(tài)[4].
圖7展示了SH0波與不同角度微裂紋作用后的切向位移時程對應的二次諧波歸一化非線性散射場(SH0-SH0模式).可以看出,二次諧波散射場隨著微裂紋取向的變化而敏感變化.具體地,對于0°微裂紋(圖7(a)),二次諧波散射場具有很好的對稱性,且只存在一個明顯的指向水平向左的主瓣,而微裂紋的取向角恰好與主瓣的對稱軸角度一致.對于20°,30°,40°,60°微裂紋(圖7(b)—7(e)),二次諧波散射場均呈現(xiàn)良好的對稱性,且微裂紋的取向角與散射場中分布于裂紋兩側的主瓣中軸線的角度幾乎一致.而對于80°微裂紋(圖7(f)),雖然二次諧波散射場依然呈現(xiàn)較好的對稱性,但是指向性相對于其他角度裂紋的散射場變弱.這是由于90°微裂紋的裂紋面平行于SH0波的質點振動方向,當裂紋面處于自由應力狀態(tài)時,即不存在法向預應力時,不會發(fā)生接觸聲非線性效應.因此,當微裂紋角度接近于90°時,接觸聲非線性效應減弱,從而導致非線性散射場的指向性會在一定程度上被破壞.即便如此,仍然可以通過散射場的對稱性對微裂紋的取向進行大致的判斷.由此可見,切向位移的二次諧波散射場可對微裂紋取向進行定量表征.本文的計算結果與文獻[30]吻合良好,證明了本文發(fā)展的數(shù)值方法的有效性和正確性.
圖7 不同微裂紋取向對應的切向位移二次諧波歸一化散射場Fig. 7 Normalized scattering fields of the 2nd harmonics for different microcrack orientations
由2.2小節(jié)的結果可知,在結構處于自由狀態(tài)時,通過SH0波與微裂紋作用的二次諧波散射場的指向性和對稱性可以很好地對微裂紋的取向進行表征.本小節(jié)繼續(xù)探討受載條件下SH0波的非線性散射場.需要注意的是,當考慮受力狀態(tài)下的非線性散射場時,整個計算過程可分為兩步進行.第一步,計算外部載荷對系統(tǒng)產生的初始應力場;第二步,在初始應力場的基礎上計算SH0波與微裂紋的非線性作用.
圖8為單軸拉伸狀態(tài)下的SH0波與不同角度微裂紋作用后的切向位移二次諧波歸一化非線性散射場.從圖中可以看出,預拉應力并不會破壞散射場的指向性和對稱性.這說明即使結構處于預拉應力狀態(tài)下,仍然可以通過二次諧波散射場來對結構中微裂紋的取向進行判斷.此外,還可以觀察到,拉應力會對非線性散射場的幅值產生明顯的影響.具體來說,二次諧波幅值隨著拉應力的增加呈下降趨勢,這是因為拉應力使微裂紋的兩個裂紋面產生初始間隙,拉應力越大,初始間隙越大,SH0波與微裂紋作用時,裂紋面周期性的接觸非線性效應減弱,從而抑制了二次諧波的產生.當拉應力足夠大時,激勵波的幅值大小無法克服此時拉應力產生的裂紋面間隙,因此不會出現(xiàn)接觸聲非線性效應,即不會產生二次諧波.
圖8 單軸拉伸狀態(tài)下不同微裂紋取向時的切向位移二次諧波歸一化散射場Fig. 8 Normalized scattering fields of the 2nd harmonics under uniaxial tension for different microcrack orientations
上述現(xiàn)象可以通過提取非線性散射場中兩個主瓣方向響應點的二次諧波幅值與拉應力之間的關系曲線來進一步說明.如圖9所示,以40°和60°微裂紋模型為例,隨著拉應力值增加,二次諧波的幅值單調遞減.隨著拉應力進一步增大,二次諧波幅值為零,即出現(xiàn)明顯的閾值效應,原因在于當SH0波的最大幅值小于該拉應力值產生的最小裂紋面間隙,將不會產生接觸聲非線性效應.
(a) 40°微裂紋 (b) 60°微裂紋(a) For the 40° microcrack (b) For the 60° microcrack圖9 主瓣方向的二次諧波幅值與拉應力的關系曲線Fig. 9 Relationships between the 2nd harmonic amplitude and the tensile stress in the direction of the main lobe
圖10為單軸壓縮狀態(tài)下的SH0波與不同角度微裂紋作用的切向位移二次諧波歸一化非線性散射場.可以觀察到,和單軸拉伸情況一樣,壓應力也不會破壞非線性散射場的指向性和對稱性,仍然可以通過散射場的對稱性來對微裂紋的取向進行判斷.此外,壓應力會對二次諧波的幅值產生顯著的影響,但是影響規(guī)律和拉應力有所不同.
圖10 單軸壓縮狀態(tài)下不同微裂紋取向對應的切向位移二次諧波歸一化散射場Fig. 10 Normalized scattering fields of the 2nd harmonics under uniaxial compression for different microcrack orientations
如圖11所示,仍然以40°和60°微裂紋為例,分別繪制非線性散射場中兩個主瓣方向響應點的二次諧波幅值與預壓應力的關系曲線.從圖中可以看到,當壓應力值增長處于較小范圍內時,二次諧波幅值隨壓應力增大而增大.隨著壓應力繼續(xù)增大,二次諧波幅值開始出現(xiàn)下降的趨勢,這是因為壓應力增大使得在裂紋面產生的接觸應力增大,因而當SH0波與微裂紋作用時,接觸聲非線性效應逐漸減弱,從而使二次諧波滋生受到抑制.當壓應力值進一步增大,SH0波的應力分量無法克服接觸應力,此時不再產生接觸聲非線性效應,二次諧波幅值變?yōu)榱悖傊?二次諧波幅值隨著壓應力的增加呈先增加后減小最后保持為零的變化趨勢,這一現(xiàn)象與二次諧波幅值隨拉應力的變化趨勢明顯不同.
(a) 40°微裂紋 (b) 60°微裂紋(a) For the 40° microcrack (b) For the 60° microcrack圖11 主瓣方向的二次諧波幅值與壓應力的關系曲線Fig. 11 Relationships between the 2nd harmonic amplitude and the compressive stress in the direction of the main lobe
本文針對SH0波與微裂紋作用的非線性散射場問題,發(fā)展了一種數(shù)值方法來處理此類問題.該方法結合了處理有限元/譜單元耦合界面的mortar方法和計算接觸摩擦問題的雙勢理論.利用mortar方法在處理非協(xié)調區(qū)域問題上的優(yōu)勢,分別采用有限單元和譜單元來離散內部的裂紋區(qū)域和外部的均勻區(qū)域,從而可以充分利用有限單元離散復雜結構的優(yōu)勢和譜單元計算效率高的優(yōu)點.在求解不同角度裂紋模型的散射場時,只需要將外部區(qū)域固定,對裂紋區(qū)域進行整體旋轉,該策略可以極大地簡化裂紋建模流程,提升建模效率,同時還保證了計算不同裂紋模型時網格的一致性.提出了一種預處理半顯式算法,該算法囊括了包含附加界面力的局部接觸力的隱式求解和全局位移的顯式求解,其可以在保證計算效率的同時不犧牲數(shù)值精度.本文的主要結論如下:
1) SH0波與微裂紋作用后的二次諧波散射場可以用于反演裂紋的取向.單軸預應力不會改變二次諧波散射場的指向性.
2) 非線性散射場中二次諧波的幅值會隨著預拉應力或預壓應力的變化而敏感變化.當結構受拉時,二次諧波幅值隨拉應力增大而單調遞減;當結構受壓時,二次諧波幅值隨壓應力增大而先增加后減小.無論何種應力狀態(tài),當應力過大時,均會出現(xiàn)閾值現(xiàn)象,超過該閾值后將不再產生非線性效應.