陳致名,余藝萌
(中國飛行試驗(yàn)研究院,陜西 西安 710089)
Gilbert[1]在20世紀(jì)70年代提出了矩張量這一物理概念,可以通過矩張量理論相關(guān)公式推導(dǎo)出地震時(shí)地下天然裂紋開裂后輻射出壓縮彈性波的信息。矩張量不僅可以描述天然地震的開裂模型,還可以在地下爆炸等領(lǐng)域描述爆炸點(diǎn)源體積的快速膨脹及由于傳播介質(zhì)相變導(dǎo)致的爆炸點(diǎn)源體積的快速坍塌[2]。Ohtsu[3]在20世紀(jì)90年代嘗試將矩張量反演理論與聲發(fā)射試驗(yàn)相結(jié)合,利用傳感器數(shù)據(jù)獲得裂紋開裂機(jī)制,但在聲發(fā)射試驗(yàn)研究中,對(duì)裂紋開裂機(jī)制的研究遇到了問題。不過隨著相關(guān)技術(shù)的進(jìn)步,傳感器精度有了明顯提高,將聲發(fā)射試驗(yàn)和矩張量反演理論相結(jié)合的相關(guān)研究也重回研究學(xué)者的視線內(nèi)[4]。Aki[5]等研究了從震源開裂形成彈性壓縮波到傳感器接收波動(dòng)數(shù)據(jù)并分析處理、反演計(jì)算的整個(gè)過程。
目前矩張量源的研究也引起了相關(guān)學(xué)者的興趣[6],在用矩張量等效體現(xiàn)震源的動(dòng)力學(xué)特性之后,僅需遠(yuǎn)場假設(shè),整體位移就可以由遠(yuǎn)場位移替代,遠(yuǎn)場位移和矩張量之間保持線性關(guān)系。地震觀測(cè)中,采用初至波幅值進(jìn)行反演計(jì)算求得震源開裂機(jī)制的等效矩張量時(shí),要求有數(shù)個(gè)位置合理的臺(tái)站可以提供臺(tái)站數(shù)據(jù),受條件因素影響,無法布置合理位置的臺(tái)站時(shí),需要從位置不佳的臺(tái)站獲得數(shù)據(jù)[7]。利用聲發(fā)射技術(shù)模擬巖石裂紋開裂時(shí),采用矩張量反演理論不僅可以得到巖石內(nèi)部裂紋面的時(shí)空演化機(jī)制,而且可以詳細(xì)描述裂紋面開裂信息,從而掌握不同裂紋開裂時(shí)裂紋之間互相干擾、互相影響的規(guī)律及單個(gè)裂紋擴(kuò)展或裂紋之間互相穿透的機(jī)制[8]。
近年來矩張量反演理論更多的應(yīng)用在水力壓裂、巖體工程、地下采礦等許多領(lǐng)域,并且在以上介紹的不同場景下,矩張量反演理論與不同實(shí)際情況相結(jié)合,有了更長足的發(fā)展[9-11]。在實(shí)際工程中,常在裂紋附近位置均勻布置傳感器,對(duì)矩張量反演計(jì)算結(jié)果的精度也無法定量判斷;而在飛機(jī)結(jié)構(gòu)中,結(jié)構(gòu)間存在互相干涉且局部可能存在大量鉚釘,由于復(fù)材與金屬結(jié)構(gòu)的共同存在,材料參數(shù)也極為復(fù)雜,利用均勻布置的傳感器矩張量反演結(jié)果精度未必可以得到保證。目前并沒有成熟的方法用來判斷復(fù)雜結(jié)構(gòu)下最優(yōu)傳感器位置,本文將有限元計(jì)算與矩張量反演相結(jié)合,提出了一種有效判斷傳感器最優(yōu)布置位置及矩張量反演精度的方法。
矩張量的物理意義是,一個(gè)微裂紋開裂的動(dòng)力學(xué)特性可以等效為在一個(gè)微體積元上加載力矩。而這個(gè)受載的微體積元其上分布的等效載荷依方向可分為9個(gè)力矩,將9個(gè)力矩按順序組合成矩陣,該矩陣即為矩張量,形式如下:
(1)
式中,右側(cè)矩陣內(nèi)的分量mij為對(duì)應(yīng)震源的等效矩張量分量。當(dāng)i與j相同時(shí),mij代表法向力矩;當(dāng)i與j不同時(shí),mij代表切向力矩。
裂紋開裂時(shí)的動(dòng)力學(xué)狀態(tài)與矩張量元素之間互相等效,公式如下:
mpq=λvk[uk(ξ,τ)]δpq+μ(vp[uq(ξ,τ)]+vq[up(ξ,τ)])
(2)
式中:mpq為矩張量元素(不同方向的等效力矩);λ、μ為拉梅常數(shù)(可由E、σ導(dǎo)出);ui為裂紋面的法向向量;vi為裂紋面的位移向量;δpq為狄拉克雷函數(shù)。
基于震源的波動(dòng)論和點(diǎn)源假設(shè),引入格林函數(shù),把矩張量與傳感器收集到的數(shù)據(jù)聯(lián)系起來。傳感器k接收到的波形位移振幅uk應(yīng)為:
uk(x,t)=Gki,j(x,t;ξ,t′)Mij(ξ,t′)
(3)
式中:*為卷積符號(hào);Gki(x,t;ξ,t′)為彈性動(dòng)力學(xué)格林函數(shù),是由單位脈沖集中力引起的位移場,即震源(ξ,t′)和監(jiān)測(cè)端(x,t)之間介質(zhì)的脈沖響應(yīng);Mij(ξ,t′)為3×3大小的矩張量矩陣。
公式(3)中Green函數(shù)的解為彈性波位移場,可以依據(jù)監(jiān)測(cè)端與裂紋中心距離的不同,將裂紋開裂輻射后的彈性波位移場的三個(gè)分量分別表示為近場項(xiàng)、中場項(xiàng)以及遠(yuǎn)場項(xiàng)。當(dāng)傳感器的位置與裂紋較遠(yuǎn)時(shí),近場項(xiàng)與中場項(xiàng)相比于遠(yuǎn)場項(xiàng)足夠小,可以將這兩項(xiàng)彈性波位移場的分量忽略,這個(gè)距離條件也被成為遠(yuǎn)場條件。
當(dāng)監(jiān)測(cè)端與裂紋之間的距離符合遠(yuǎn)場條件時(shí),矩張量反演公式可以簡化為下式:
(4)
式(4)為三維空間下的矩張量反演公式。由于本文的研究內(nèi)容范圍為二維空間,因此需將此方程進(jìn)一步簡化:
(5)
在有限元計(jì)算過程中,需要構(gòu)造一個(gè)簡單的裂紋模型,篩選出符合遠(yuǎn)場條件的節(jié)點(diǎn)編號(hào),再求出裂紋開裂后以上節(jié)點(diǎn)的初至波幅值,將幅值帶入矩張量反演計(jì)算,通過對(duì)比得到最優(yōu)的傳感器布置位置。
本文選用物理模型為半徑400 mm的圓,將傳感器布置在半徑300 mm的圓上,點(diǎn)源設(shè)置在正中心,傳感器與點(diǎn)源距離大于波長的8倍,滿足遠(yuǎn)場條件(由于飛機(jī)上被監(jiān)測(cè)結(jié)構(gòu)之間形狀差別很大,因此模型的形狀并無具體要求,只需滿足傳感器的布置條件即可,模型具體形狀依據(jù)劃分網(wǎng)格的方式而定)。
圖1 物理模型
圖1為一種物理模型示意圖。
由于裂紋開裂過程為動(dòng)態(tài)受載,整個(gè)加載過程中力的變化類似于階躍載荷。因此在有限元建模時(shí)對(duì)裂紋面施加強(qiáng)制位移來模擬裂紋的斷裂過程,經(jīng)歸一化處理,強(qiáng)制位移的最大值設(shè)為1,加載點(diǎn)的強(qiáng)制位移與時(shí)間的函數(shù)關(guān)系如式(6)所示:
(6)
式中:Tr為半周期長,時(shí)間到達(dá)Tr時(shí)位移的幅值達(dá)到峰值。強(qiáng)制位移-時(shí)間函數(shù)圖如圖2所示。
圖2 強(qiáng)制位移-時(shí)間函數(shù)圖
首先將有限元計(jì)算結(jié)果整理成極坐標(biāo)的形式與理論結(jié)果比較,其后從所有符合條件的節(jié)點(diǎn)中任取三個(gè)節(jié)點(diǎn)反演計(jì)算出點(diǎn)源矩張量,將其與設(shè)定值進(jìn)行對(duì)比,利用蒙特卡洛方法進(jìn)行大量重復(fù)試驗(yàn),從眾多樣本中篩選出一系列使計(jì)算結(jié)果誤差最小的傳感器布置組合,將每個(gè)組合連線成三角形,并畫在一張圖上,就可以最優(yōu)的傳感器布置位置的分布方式,同時(shí)把得到最優(yōu)傳感器組合的樣本數(shù)與總樣本數(shù)進(jìn)行比較,比值的高低可視為反演精度的高低。
得到反演矩張量后,需要對(duì)其精度進(jìn)行判斷,本文給出一種判斷方法:首先通過矩張量反演計(jì)算從有限元計(jì)算結(jié)果中得到計(jì)算矩張量,將其進(jìn)行歸一化處理,分別計(jì)算該矩陣與有限元建模時(shí)的加載等效矩張量的Frobenius范數(shù),將兩范數(shù)比較得到誤差,同時(shí)還需比較計(jì)算矩張量矩陣與加載等效矩張量矩陣的主元素大小差別,最后只有當(dāng)范數(shù)誤差與主元素誤差同時(shí)小于一定程度時(shí),才可以將該組合視為最優(yōu)組合。
分別采用四節(jié)點(diǎn)四邊形單元和三節(jié)點(diǎn)三角形單元兩種有限元網(wǎng)格劃分形式分析矩張量反演問題,并比較了兩種簡單裂紋模型的精度區(qū)別。
模型材料參數(shù)設(shè)定為:E=5.4 GPa,σ=0.2,ρ=3000 kg/m3。
裂紋的構(gòu)造方式為在同一位置放置兩個(gè)節(jié)點(diǎn),兩點(diǎn)之間不作約束,同時(shí)對(duì)兩個(gè)節(jié)點(diǎn)施加方向不同的強(qiáng)制位移。
2.1.1 模型介紹
(7)
2.1.2 計(jì)算驗(yàn)證
首先整理所需節(jié)點(diǎn)編號(hào),在NASTRAN中導(dǎo)出所有節(jié)點(diǎn)的位移文件后,提取初至波峰值,同時(shí)利用矩張量反演公式,可以求得任意節(jié)點(diǎn)的理論初至波峰值,將兩者以極坐標(biāo)的形式進(jìn)行比較并歸一化。
由圖3和圖4可知,對(duì)于拉伸裂紋來說,三角形網(wǎng)格的計(jì)算結(jié)果與理論結(jié)果吻合效果較好,整體無較大計(jì)算誤差。
圖3 四邊形網(wǎng)格對(duì)比圖 圖4 三角形網(wǎng)格對(duì)比圖
2.2.1 模型介紹
(8)
2.2.2 計(jì)算驗(yàn)證
計(jì)算初至波峰值與理論初至波峰值比較圖如圖5和圖6所示。
圖5 三角形網(wǎng)格對(duì)比圖 圖6 四邊形網(wǎng)格對(duì)比圖
由圖5和圖6可知,對(duì)于剪切裂紋來說,三角形網(wǎng)格的計(jì)算結(jié)果與理論結(jié)果吻合效果較好,整體無較大計(jì)算誤差。
矩張量反演問題與常規(guī)裂紋計(jì)算問題不同,矩張量反演問題關(guān)注不同方向幅值的相對(duì)大小,而不是有限元計(jì)算結(jié)果與真實(shí)結(jié)果之間的絕對(duì)大小。三角形單元作為常應(yīng)變單元計(jì)算精度較低,但是其在單元的不同邊上計(jì)算位移與真實(shí)位移之間的誤差近似,誤差的空間分布更平均,因此用三角形單元有可能會(huì)得到更好的矩張量計(jì)算結(jié)果。
兩種裂紋的理論幅值與計(jì)算幅值的擬合效果都很好。而由于拉伸裂紋的幅值形狀為雙頁,剪切裂紋的幅值形狀為四頁,雙頁形狀變化更便于觀察,因此在分析材料屬性及材料缺陷產(chǎn)生的影響時(shí),選擇拉伸裂紋模型。
本章分析材料參數(shù)對(duì)反演精度產(chǎn)生的影響,將不同材料參數(shù)導(dǎo)入模型進(jìn)行矩張量反演計(jì)算,分析矩張量反演方法應(yīng)用于各種不同材料時(shí)的精度變化。
分別對(duì)不同泊松比模型分析計(jì)算,比較各個(gè)模型符合反演精度的最優(yōu)組合占比。
表1 反演精度隨泊松比變化
泊松比為0.2時(shí)最優(yōu)組合分布如圖7所示(將每個(gè)最優(yōu)組合用三角形表示)。
圖7 泊松比為0.2時(shí)最優(yōu)組合分布圖
可得最優(yōu)組合的位置分布并無明顯規(guī)律,理想拉伸裂紋模型的最優(yōu)傳感器選擇受較低限制,最優(yōu)組合占比低于100%的原因是篩除掉了傳感器位置較近的傳感器組合,且將最優(yōu)組合的誤差條件設(shè)置得較高。
由于泊松比影響較小,因此有限元計(jì)算時(shí),不用過多考慮泊松比的變化對(duì)反演精度的影響。
分別對(duì)不同彈性模量的模型分析計(jì)算,比較各個(gè)模型符合反演精度的最優(yōu)組合占比。
表2 反演精度隨彈性模量變化
隨著彈性模量的變化,反演精度變化較小。當(dāng)彈性模量變大時(shí),最優(yōu)組合數(shù)量占比略為減小。
選用兩種復(fù)合材料進(jìn)行分析計(jì)算,兩種復(fù)合材料的材料參數(shù)見表3,幅值比較圖如圖8和圖9所示。
表3 兩種復(fù)合材料的材料參數(shù)
由圖8和圖9可知,兩種復(fù)合材料的計(jì)算與理論幅值曲線圖相差不大,其最優(yōu)組合占比見表4。
圖8 復(fù)材1對(duì)比圖 圖9 復(fù)材2對(duì)比圖
表4 不同各向異性材料的最優(yōu)組合占比
復(fù)材1的最優(yōu)組合分布如圖10所示。
圖10 復(fù)材1的最優(yōu)組合分布
可得以上兩種復(fù)合材料與理想模型相比反演精度差別較小,且最優(yōu)傳感器位置分布規(guī)律并未因材料為復(fù)合材料而發(fā)生明顯變化,即對(duì)于常見的復(fù)合材料矩張量反演計(jì)算時(shí),不需要特意調(diào)整傳感器布置位置。
1)有限元方法配合蒙特卡羅方法可以有效輔助矩張量反演問題的計(jì)算,篩選總結(jié)出最優(yōu)傳感器位置;
2)有限元模型網(wǎng)格劃分形式對(duì)矩張量反演計(jì)算結(jié)果的影響較大,單元精度的空間均勻性比絕對(duì)大小更重要;
3)簡單模型的最優(yōu)傳感器位置幾乎不受材料參數(shù)影響,對(duì)于飛機(jī)的復(fù)雜結(jié)構(gòu)(如多條裂紋同時(shí)開裂、局部存在多個(gè)鉚釘、結(jié)構(gòu)存在孔隙及孔洞等),可通過本文提出的方法有效判斷其傳感器的最優(yōu)位置。