基于有限元研究傳感器位置對(duì)矩張量反演方法精度影響

陳致名，余藝萌

(中國(guó)飛行試驗(yàn)研究院，陜西 西安 710089)

0 引言

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)場(chǎng)假設(shè)，整體位移就可以由遠(yuǎn)場(chǎng)位移替代，遠(yuǎn)場(chǎng)位移和矩張量之間保持線性關(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)域，并且在以上介紹的不同場(chǎng)景下，矩張量反演理論與不同實(shí)際情況相結(jié)合，有了更長(zhǎng)足的發(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)布置位置及矩張量反演精度的方法。

1 研究方法

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

矩張量的物理意義是，一個(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ù)，是由單位脈沖集中力引起的位移場(chǎng)，即震源(ξ,t′)和監(jiān)測(cè)端(x，t)之間介質(zhì)的脈沖響應(yīng)；Mij(ξ,t′)為3×3大小的矩張量矩陣。

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

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

(4)

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

(5)

1.2 有限元計(jì)算內(nèi)容

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

本文選用物理模型為半徑400 mm的圓，將傳感器布置在半徑300 mm的圓上，點(diǎn)源設(shè)置在正中心，傳感器與點(diǎn)源距離大于波長(zhǎng)的8倍，滿足遠(yuǎn)場(chǎng)條件(由于飛機(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為半周期長(zhǎng)，時(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)組合。

2 簡(jiǎn)單模型計(jì)算結(jié)果

分別采用四節(jié)點(diǎn)四邊形單元和三節(jié)點(diǎn)三角形單元兩種有限元網(wǎng)格劃分形式分析矩張量反演問題，并比較了兩種簡(jiǎn)單裂紋模型的精度區(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 拉伸裂紋模型

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 剪切裂紋模型

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ì)算誤差。

2.3 總結(jié)

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

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

3 材料參數(shù)的影響

本章分析材料參數(shù)對(duì)反演精度產(chǎn)生的影響，將不同材料參數(shù)導(dǎo)入模型進(jìn)行矩張量反演計(jì)算，分析矩張量反演方法應(yīng)用于各種不同材料時(shí)的精度變化。

3.1 泊松比

分別對(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ì)反演精度的影響。

3.2 彈性模量

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

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

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

3.3 復(fù)合材料的適用性

選用兩種復(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)整傳感器布置位置。

4 結(jié)論

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

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

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