陳致名，余藝萌

(中國飛行試驗研究院，陜西 西安 710089)

0 引言

Gilbert[1]在20世紀70年代提出了矩張量這一物理概念，可以通過矩張量理論相關公式推導出地震時地下天然裂紋開裂后輻射出壓縮彈性波的信息。矩張量不僅可以描述天然地震的開裂模型，還可以在地下爆炸等領域描述爆炸點源體積的快速膨脹及由于傳播介質相變導致的爆炸點源體積的快速坍塌[2]。Ohtsu[3]在20世紀90年代嘗試將矩張量反演理論與聲發(fā)射試驗相結合，利用傳感器數據獲得裂紋開裂機制，但在聲發(fā)射試驗研究中，對裂紋開裂機制的研究遇到了問題。不過隨著相關技術的進步，傳感器精度有了明顯提高，將聲發(fā)射試驗和矩張量反演理論相結合的相關研究也重回研究學者的視線內[4]。Aki[5]等研究了從震源開裂形成彈性壓縮波到傳感器接收波動數據并分析處理、反演計算的整個過程。

目前矩張量源的研究也引起了相關學者的興趣[6]，在用矩張量等效體現震源的動力學特性之后，僅需遠場假設，整體位移就可以由遠場位移替代，遠場位移和矩張量之間保持線性關系。地震觀測中，采用初至波幅值進行反演計算求得震源開裂機制的等效矩張量時，要求有數個位置合理的臺站可以提供臺站數據，受條件因素影響，無法布置合理位置的臺站時，需要從位置不佳的臺站獲得數據[7]。利用聲發(fā)射技術模擬巖石裂紋開裂時，采用矩張量反演理論不僅可以得到巖石內部裂紋面的時空演化機制，而且可以詳細描述裂紋面開裂信息，從而掌握不同裂紋開裂時裂紋之間互相干擾、互相影響的規(guī)律及單個裂紋擴展或裂紋之間互相穿透的機制[8]。

近年來矩張量反演理論更多的應用在水力壓裂、巖體工程、地下采礦等許多領域，并且在以上介紹的不同場景下，矩張量反演理論與不同實際情況相結合，有了更長足的發(fā)展[9-11]。在實際工程中，常在裂紋附近位置均勻布置傳感器，對矩張量反演計算結果的精度也無法定量判斷；而在飛機結構中，結構間存在互相干涉且局部可能存在大量鉚釘，由于復材與金屬結構的共同存在，材料參數也極為復雜，利用均勻布置的傳感器矩張量反演結果精度未必可以得到保證。目前并沒有成熟的方法用來判斷復雜結構下最優(yōu)傳感器位置，本文將有限元計算與矩張量反演相結合，提出了一種有效判斷傳感器最優(yōu)布置位置及矩張量反演精度的方法。

1 研究方法

1.1 矩張量理論和矩張量反演理論

矩張量的物理意義是，一個微裂紋開裂的動力學特性可以等效為在一個微體積元上加載力矩。而這個受載的微體積元其上分布的等效載荷依方向可分為9個力矩，將9個力矩按順序組合成矩陣，該矩陣即為矩張量，形式如下：

(1)

式中，右側矩陣內的分量mij為對應震源的等效矩張量分量。當i與j相同時，mij代表法向力矩；當i與j不同時，mij代表切向力矩。

裂紋開裂時的動力學狀態(tài)與矩張量元素之間互相等效，公式如下：

mpq=λvk[uk(ξ,τ)]δpq+μ(vp[uq(ξ,τ)]+vq[up(ξ,τ)])

(2)

式中：mpq為矩張量元素(不同方向的等效力矩)；λ、μ為拉梅常數(可由E、σ導出)；ui為裂紋面的法向向量；vi為裂紋面的位移向量；δpq為狄拉克雷函數。

基于震源的波動論和點源假設，引入格林函數，把矩張量與傳感器收集到的數據聯(lián)系起來。傳感器k接收到的波形位移振幅uk應為：

uk(x,t)=Gki，j(x,t;ξ,t′)Mij(ξ,t′)

(3)

式中：*為卷積符號；Gki(x,t;ξ,t′)為彈性動力學格林函數，是由單位脈沖集中力引起的位移場，即震源(ξ,t′)和監(jiān)測端(x，t)之間介質的脈沖響應；Mij(ξ,t′)為3×3大小的矩張量矩陣。

公式(3)中Green函數的解為彈性波位移場，可以依據監(jiān)測端與裂紋中心距離的不同，將裂紋開裂輻射后的彈性波位移場的三個分量分別表示為近場項、中場項以及遠場項。當傳感器的位置與裂紋較遠時，近場項與中場項相比于遠場項足夠小，可以將這兩項彈性波位移場的分量忽略，這個距離條件也被成為遠場條件。

當監(jiān)測端與裂紋之間的距離符合遠場條件時，矩張量反演公式可以簡化為下式：

(4)

式(4)為三維空間下的矩張量反演公式。由于本文的研究內容范圍為二維空間，因此需將此方程進一步簡化：

(5)

1.2 有限元計算內容

在有限元計算過程中，需要構造一個簡單的裂紋模型，篩選出符合遠場條件的節(jié)點編號，再求出裂紋開裂后以上節(jié)點的初至波幅值，將幅值帶入矩張量反演計算，通過對比得到最優(yōu)的傳感器布置位置。

本文選用物理模型為半徑400 mm的圓，將傳感器布置在半徑300 mm的圓上，點源設置在正中心，傳感器與點源距離大于波長的8倍，滿足遠場條件(由于飛機上被監(jiān)測結構之間形狀差別很大，因此模型的形狀并無具體要求，只需滿足傳感器的布置條件即可，模型具體形狀依據劃分網格的方式而定)。

圖1 物理模型

圖1為一種物理模型示意圖。

由于裂紋開裂過程為動態(tài)受載，整個加載過程中力的變化類似于階躍載荷。因此在有限元建模時對裂紋面施加強制位移來模擬裂紋的斷裂過程，經歸一化處理，強制位移的最大值設為1，加載點的強制位移與時間的函數關系如式(6)所示：

(6)

式中：Tr為半周期長，時間到達Tr時位移的幅值達到峰值。強制位移-時間函數圖如圖2所示。

圖2 強制位移-時間函數圖

首先將有限元計算結果整理成極坐標的形式與理論結果比較，其后從所有符合條件的節(jié)點中任取三個節(jié)點反演計算出點源矩張量，將其與設定值進行對比，利用蒙特卡洛方法進行大量重復試驗，從眾多樣本中篩選出一系列使計算結果誤差最小的傳感器布置組合，將每個組合連線成三角形，并畫在一張圖上，就可以最優(yōu)的傳感器布置位置的分布方式，同時把得到最優(yōu)傳感器組合的樣本數與總樣本數進行比較，比值的高低可視為反演精度的高低。

得到反演矩張量后，需要對其精度進行判斷，本文給出一種判斷方法：首先通過矩張量反演計算從有限元計算結果中得到計算矩張量，將其進行歸一化處理，分別計算該矩陣與有限元建模時的加載等效矩張量的Frobenius范數，將兩范數比較得到誤差，同時還需比較計算矩張量矩陣與加載等效矩張量矩陣的主元素大小差別，最后只有當范數誤差與主元素誤差同時小于一定程度時，才可以將該組合視為最優(yōu)組合。

2 簡單模型計算結果

分別采用四節(jié)點四邊形單元和三節(jié)點三角形單元兩種有限元網格劃分形式分析矩張量反演問題，并比較了兩種簡單裂紋模型的精度區(qū)別。

模型材料參數設定為：E=5.4 GPa，σ=0.2，ρ=3000 kg/m3。

裂紋的構造方式為在同一位置放置兩個節(jié)點，兩點之間不作約束，同時對兩個節(jié)點施加方向不同的強制位移。

2.1 拉伸裂紋模型

2.1.1 模型介紹

(7)

2.1.2 計算驗證

首先整理所需節(jié)點編號，在NASTRAN中導出所有節(jié)點的位移文件后，提取初至波峰值，同時利用矩張量反演公式，可以求得任意節(jié)點的理論初至波峰值，將兩者以極坐標的形式進行比較并歸一化。

由圖3和圖4可知，對于拉伸裂紋來說，三角形網格的計算結果與理論結果吻合效果較好，整體無較大計算誤差。

圖3 四邊形網格對比圖 圖4 三角形網格對比圖

2.2 剪切裂紋模型

2.2.1 模型介紹

(8)

2.2.2 計算驗證

計算初至波峰值與理論初至波峰值比較圖如圖5和圖6所示。

圖5 三角形網格對比圖 圖6 四邊形網格對比圖

由圖5和圖6可知，對于剪切裂紋來說，三角形網格的計算結果與理論結果吻合效果較好，整體無較大計算誤差。

2.3 總結

矩張量反演問題與常規(guī)裂紋計算問題不同，矩張量反演問題關注不同方向幅值的相對大小，而不是有限元計算結果與真實結果之間的絕對大小。三角形單元作為常應變單元計算精度較低，但是其在單元的不同邊上計算位移與真實位移之間的誤差近似，誤差的空間分布更平均，因此用三角形單元有可能會得到更好的矩張量計算結果。

兩種裂紋的理論幅值與計算幅值的擬合效果都很好。而由于拉伸裂紋的幅值形狀為雙頁，剪切裂紋的幅值形狀為四頁，雙頁形狀變化更便于觀察，因此在分析材料屬性及材料缺陷產生的影響時，選擇拉伸裂紋模型。

3 材料參數的影響

本章分析材料參數對反演精度產生的影響，將不同材料參數導入模型進行矩張量反演計算，分析矩張量反演方法應用于各種不同材料時的精度變化。

3.1 泊松比

分別對不同泊松比模型分析計算，比較各個模型符合反演精度的最優(yōu)組合占比。

表1 反演精度隨泊松比變化

泊松比為0.2時最優(yōu)組合分布如圖7所示(將每個最優(yōu)組合用三角形表示)。

圖7 泊松比為0.2時最優(yōu)組合分布圖

可得最優(yōu)組合的位置分布并無明顯規(guī)律，理想拉伸裂紋模型的最優(yōu)傳感器選擇受較低限制，最優(yōu)組合占比低于100%的原因是篩除掉了傳感器位置較近的傳感器組合，且將最優(yōu)組合的誤差條件設置得較高。

由于泊松比影響較小，因此有限元計算時，不用過多考慮泊松比的變化對反演精度的影響。

3.2 彈性模量

分別對不同彈性模量的模型分析計算，比較各個模型符合反演精度的最優(yōu)組合占比。

表2 反演精度隨彈性模量變化

隨著彈性模量的變化，反演精度變化較小。當彈性模量變大時，最優(yōu)組合數量占比略為減小。

3.3 復合材料的適用性

選用兩種復合材料進行分析計算，兩種復合材料的材料參數見表3，幅值比較圖如圖8和圖9所示。

表3 兩種復合材料的材料參數

由圖8和圖9可知，兩種復合材料的計算與理論幅值曲線圖相差不大，其最優(yōu)組合占比見表4。

圖8 復材1對比圖 圖9 復材2對比圖

表4 不同各向異性材料的最優(yōu)組合占比

復材1的最優(yōu)組合分布如圖10所示。

圖10 復材1的最優(yōu)組合分布

可得以上兩種復合材料與理想模型相比反演精度差別較小，且最優(yōu)傳感器位置分布規(guī)律并未因材料為復合材料而發(fā)生明顯變化，即對于常見的復合材料矩張量反演計算時，不需要特意調整傳感器布置位置。

4 結論

1)有限元方法配合蒙特卡羅方法可以有效輔助矩張量反演問題的計算，篩選總結出最優(yōu)傳感器位置；

2)有限元模型網格劃分形式對矩張量反演計算結果的影響較大，單元精度的空間均勻性比絕對大小更重要；

3)簡單模型的最優(yōu)傳感器位置幾乎不受材料參數影響，對于飛機的復雜結構(如多條裂紋同時開裂、局部存在多個鉚釘、結構存在孔隙及孔洞等)，可通過本文提出的方法有效判斷其傳感器的最優(yōu)位置。