魯光銀， 張志勇.2， 朱自強， 付渝渝

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長沙 410083; 2.河南省交通規(guī)劃設(shè)計研究陳院股份有限公司，鄭州 450000)

0 引言

近年來，隨著各級建設(shè)工程要求地提高，混凝土成為各種工程建筑最常用的一種結(jié)構(gòu)材料[1]。在混凝土施工過程中，由于混凝土材質(zhì)、外界因素的影響，在內(nèi)部結(jié)構(gòu)中往往會有缺陷存在[2]，但混凝土的質(zhì)量對整個結(jié)構(gòu)工程至關(guān)重要，關(guān)乎到工程的使用安全性和使用期限[3]。因此，在混凝土澆筑完成之后，需要對混凝土的質(zhì)量進(jìn)行相應(yīng)的檢測，確保其安全性。同時，在對混凝土檢測時，又要盡量保證混凝土原有的質(zhì)量盡可能地不受到破壞[4]。因此，常用無損檢測進(jìn)行混凝土質(zhì)量的判斷。而超聲波具有很強的穿透能力，方向性強，檢測設(shè)備簡單，儀器操作方便，節(jié)省人力和檢測成本，分辨率較高，因此在混凝土無損檢測領(lǐng)域應(yīng)用最為廣泛[5]。

為了更確切地研究超聲在混凝土中的傳播規(guī)律，不同的數(shù)值計算方法逐漸被引用到超聲數(shù)值模擬中來。模擬超聲波傳播常用的數(shù)值方法有：有限元法、彈性動力學(xué)有限積分技術(shù)、有限差分法等[6]。

筆者主要介紹有限差分法在混凝土正演中的應(yīng)用進(jìn)展。有限差分法在彈性波數(shù)值模擬方面應(yīng)用最為廣泛，因此逐漸地被應(yīng)用到超聲波的數(shù)值模擬中。Day等[7]使用Pade近似法進(jìn)行了二維時間域的數(shù)值模擬；Emmerich等[8]提出了廣義標(biāo)準(zhǔn)線性固體，推導(dǎo)了粘彈性模量的有理式;Robertsson等[9]提出了速度一應(yīng)力有限交錯網(wǎng)格法，時間上二階精度，空間上四階精度，用這種方法研究了二維空間中橫波和縱波的傳播性質(zhì)與規(guī)律，并從二維推廣到三維，取得了較好的效果；H.Yamawaki等[10]在常規(guī)的有限差分交錯網(wǎng)格的基礎(chǔ)上，提出了改進(jìn)交錯網(wǎng)格節(jié)點的計算方案，使得有限差分法更能適應(yīng)于各項異性介質(zhì)；魏東等[11]應(yīng)用時域有限差分法，研究了彈性固體中的二維聲場特性，采用PML(完全匹配層)吸正收邊界條件，解決了有限區(qū)域的計算問題；張俊等[12]基于有限差分法，設(shè)置建立了含有多個探頭的相控陣超聲檢測模型；密士文[13]在混凝土超聲數(shù)值模擬波動方程中加入粘彈性介質(zhì)的拉梅常數(shù)項。但其只是推導(dǎo)出速度—應(yīng)力公式，未對超聲正演模擬效果做出相應(yīng)的研究。

相對于彈性介質(zhì)，其受到應(yīng)力和產(chǎn)生應(yīng)變是瞬發(fā)的，且滿足Hook定律中應(yīng)變與應(yīng)力成線性關(guān)系[14]。并且，當(dāng)失去外部應(yīng)力后，彈性體會迅速恢復(fù)到原來的狀態(tài)。而對于粘彈性介質(zhì)，其物理性質(zhì)介于彈性與粘性之間。在受到外力作用時，粘彈性介質(zhì)會永久發(fā)生形變，且在外力消失后，不會恢復(fù)到原來的狀態(tài)。這種性質(zhì)不僅跟粘彈性介質(zhì)自身因素有關(guān)，還與壓力大小、溫度等外界因素有關(guān)。粘彈性介質(zhì)對混凝土中能量的衰減機理比較復(fù)雜。

混凝土作為一種多項非均勻復(fù)合材料，超聲在其內(nèi)部傳播，在做數(shù)值模擬時，彈性動力方程已經(jīng)不能更確切地模擬超聲聲波的傳播方式。筆者分析總結(jié)前人工作，把有限差分法應(yīng)用到混凝土超聲波數(shù)值模擬上來。將混凝土視為粘彈性介質(zhì)，基于Kelvin粘彈性介質(zhì)[15]，詳細(xì)推導(dǎo)出了超聲在粘彈性介質(zhì)中的聲波方程，在時間二階、空間任意偶數(shù)精度下的交錯網(wǎng)格差分格式，利用數(shù)值模擬的方法來模擬混凝土中超聲波的傳播過程。

1 混凝土中超聲波傳播的數(shù)值模擬

粘彈性體常見的基本模型有：Maxwell 模型、Kelvin模型、標(biāo)準(zhǔn)線性模型和達(dá)朗貝爾模型。

目前在處理粘彈性介質(zhì)中波動方程問題主要使用兩種方法：①復(fù)模量法(粘性參數(shù)法)；②添加粘彈性介質(zhì)品質(zhì)因子Q值的方法。在全波列模擬中，通常采用的是粘性參數(shù)法，而在波場、波形分析時，則常采用Q值法[16]。筆者主要進(jìn)行聲波傳播某一時刻波場快照分析和提取某一道聲波信號進(jìn)行分析，因此，采用的是添加Q值的方法。

1.1 Kelvin粘彈性介質(zhì)

二維的Kelvin粘彈性聲波方程的一階速度—應(yīng)力方程表示形式為：

(1)

式(1)中的拉梅系數(shù)計算式為：

λ+2μ=ρv2

(2)

(3)

其中:λ和μ表示彈性介質(zhì)的拉梅系數(shù)；λ′和μ′表示粘性介質(zhì)的拉梅系數(shù)；ν表示模型介質(zhì)的縱波速度；ω為圓頻率。

在粘彈性介質(zhì)中，用Q表示縱波的品質(zhì)因子，反映模型介質(zhì)對于超聲波能量的吸收情況。式(4)中，f表示超聲波的頻率，C為超聲波的相速度，a表示超聲波振幅隨距離的衰減系數(shù)。品質(zhì)因子Q與頻率成函數(shù)關(guān)系，不同頻率的超聲波的能量衰減情況不同[16]。頻率越高，能量衰減越大，相對應(yīng)的Q值越大。但品質(zhì)因子Q的值表示的是粘彈性介質(zhì)本身的能量衰減屬性，只與粘彈性介質(zhì)物理屬性有關(guān)，與頻率的大小無關(guān)，因此，在進(jìn)行數(shù)值模擬時，通常是采用Q值固定的方法。品質(zhì)因子Q值越大，表示介質(zhì)的粘滯性越小，即介質(zhì)屬性越趨近于彈性。當(dāng)Q值為無窮大時，式(3)可表示為λ′+2μ′=0，即變成了完全彈性介質(zhì)聲波方程的表達(dá)形式。

1.2 高階精度差分粘彈性聲波方程

采用交錯網(wǎng)格法計算粘彈性聲波方程，分別取質(zhì)點在和的應(yīng)力和速度的值進(jìn)行計算。得到速度—應(yīng)力在時間一階導(dǎo)數(shù)上的二階精度差分格式：

對式(4)進(jìn)行交錯網(wǎng)絡(luò)剖分，i表示網(wǎng)格計算節(jié)點的深度方向序號，j表示節(jié)點水平方向序號，k表示時間方向序號，i+1/2、j+1/2、k+1/2分別表示相對應(yīng)的半程節(jié)點序號。分別在時間節(jié)點k和空間節(jié)點j+1/2位置處進(jìn)行計算水平速度分量vx和垂直速度分量vz。

{px(t+Δt2)=px(t-Δt2)+(λ+2μ)?vx?x+(λ′+2μ′)?vx?x

pz(t+Δt2)=pz(t-Δt2)+(λ+2μ)?vz?z+(λ′+2μ′)?vz?z

vx(t)=vx(t-Δt)+Δtρ?p?x

vz(t)=vz(t-Δt)+Δtρ?p?z

(4)

為了便于說明，在差分格式中，將水平速度分量vx設(shè)為U，垂直速度分量vz設(shè)為W；水平加速度分量?U?t=?vx?t設(shè)為U1，垂直加速度分量?W?t=?vz?t設(shè)為W1；設(shè)P、S、T分別表示應(yīng)力p、水平應(yīng)力分量px和垂直應(yīng)力分量pz。得出的時間2階，空間2N階精度的交錯網(wǎng)絡(luò)粘彈性質(zhì)聲波方程的差分格式表示為：

水平應(yīng)力分量為：

S|k+12i,j=S|k-12i,j+Δt(λ+2μ)(i,j)

1Δx∑Nn=1CNn[Uk(i,j+2n-12)-Uk(i,j-2n-12)]+

Δt(λ' + 2μ')(i,j)1Δx∑Nn=1CNn

[Uk1(i,j+2n-12)-Uk1(i,j-2n-12)]

(5)

垂直應(yīng)力分量：

T|k + 12i,j=T|k-12i,j+ Δt(λ+ 2μ)(i,j)

1Δz∑Nn=1CNn[Wk(i+2n-12,j)-Wk(i-2n-12,j)]+

Δt(λ' + 2μ')(i,j)1Δz∑Nn=1CNn

[Wk1(i + 2n-12,j)-Wk1(i-2n-12,j)]

(6)

水平速度分量：

U|ki,j+12=U|k-1i,j+12+1ρi,j+12

ΔtΔx∑Nn=1CNn(Pki,j+n-Pki,j-n+1)

(7)

垂直速度分量：

W|ki+12,j=W|k-1i+12,j+1ρi+12,j

ΔtΔz∑Nn=1CNn(Pki+n,j-Pki-n+1,j)

(8)

水平加速度分量：

U1|ki,j+12=1ρi,j+12ΔtΔx∑Nn=1CNn(Pki,j+n-Pki,j-n+1)

(9)

垂直加速度分量：

W1|ki+12,j=1ρi+12,jΔtΔz∑Nn=1CNn(Pki+n,j-Pki-n+1,j)

(10)

1.3 完全匹配層吸收邊界

圖1 穩(wěn)定系數(shù)與空間差分精度關(guān)系曲線Fig.1 Relationship between stability coefficient and spatial differential accuracy

采用完全匹配層吸收邊界(PML),在計算區(qū)域四周區(qū)域選取50層差分網(wǎng)格厚度的吸收層，用來吸收和衰減向模型之外傳播的波，假設(shè)在建立的模型區(qū)域與匹配層的交界面處產(chǎn)生最小的虛假反射[17]。

1.4 穩(wěn)定性條件

在時間2M階、空間2N階精度下，二維聲波交錯網(wǎng)格形式下的有限差分方程的穩(wěn)定性條件為：

0≤∑Mm=1{(-1)m-1(2m-1)!(Δt2S)2m·

[∑Nn=1CNn(-1)n-1]2m}≤1

(11)

差分精度選取的階數(shù)不同，穩(wěn)定性也不同，精度要求越高，穩(wěn)定性條件選取的越嚴(yán)格[18]。當(dāng)空間差分精度達(dá)到12階時，穩(wěn)定系數(shù)變化值已經(jīng)很小，再增加空間差分精度對穩(wěn)定系數(shù)的影響不大，因此，本文在做有限差分計算時，選取空間12階差分格式。

2 數(shù)值模擬分析

通過對混凝土模型進(jìn)行定量化和數(shù)字化，對混凝土中的超聲波場進(jìn)行計算，模擬聲波傳播特性。通過數(shù)值模擬可以真實有效的重復(fù)計算超聲在混凝土中的傳播過程及變化情況，分析超聲波頻率、相位、振幅及波的類型、能量等物性參數(shù)的變化，為推斷混凝土內(nèi)部結(jié)構(gòu)可能存在的缺陷提供借鑒。

2.1 混凝土骨料模型

骨料是混凝土中除水泥漿外所占體積比最大的成分，骨料形狀不規(guī)則，分布無規(guī)律。模型尺寸設(shè)置為400 mm×200 mm規(guī)格，骨料最大顆粒直徑為Amax=14 mm(圖2)。

從圖3可以看出：在均勻水泥砂漿介質(zhì)模型中，振幅衰減曲線平滑，振幅衰減值與傳播距離基本成正比關(guān)系，可以推測，均勻介質(zhì)中能量衰減應(yīng)該只與波前擴散的距離有關(guān)。在含骨料混凝土介質(zhì)中，衰減曲線不平滑，衰減值無規(guī)律，較均勻水泥砂漿介質(zhì)衰減更快，這是由于骨料分布無規(guī)律、界面不規(guī)則，超聲在骨料界面發(fā)生反射、透射及多重散射，能量衰減程度也無法確切判斷。圖4為兩種介質(zhì)模型(200,50)處質(zhì)點振動對比圖，從圖4中可以看出：兩種介質(zhì)中質(zhì)點振動速度不同，且混凝土中的質(zhì)點相位相對比較滯后，這是由于混凝土自身的粘彈性，使得介質(zhì)中質(zhì)點在震源波消失后，仍會處于振動狀態(tài)，不會立即靜止，在這一時間段內(nèi)存在剩余的應(yīng)變，不僅會造成能量的衰減，還會對信號產(chǎn)生干擾。

圖2 混凝土模型Fig.2 Concrete model

圖3 水泥漿、混凝土衰減對比Fig.3 The attenuation comparison between cement slurry and concrete

圖4 質(zhì)點振動圖Fig.4 Particle vibration map

圖5 波場快照對比圖Fig.5 The wavefield snapshot(a)t=24 μs時刻波場快照；(b)t=36 μs時刻波場快照

圖5為超聲波在含有骨料的混凝土模型中傳播不同時刻的波場快照。如圖5(a)所示，是聲波縱波、橫波的初至波波前面的位置，此刻超聲波包絡(luò)面相對比較規(guī)則連續(xù)；圖5(b)所示，聲波傳播范圍內(nèi)，由于受不規(guī)則骨料的影響，超聲在骨料界面發(fā)生多重散射、發(fā)射等現(xiàn)象，造成波形包絡(luò)面不連續(xù)，且包絡(luò)面內(nèi)波形出現(xiàn)紊亂。為了更加清晰地分析超聲傳播過程中能量的衰減情況，抽取了混凝土模型橫向中心位置(200,50)處的單道波形，如圖6所示：超聲信號整體上振幅逐漸減小，與波前擴散的衰減有關(guān)；但振幅衰減不規(guī)則變化，這是由于在混凝土介質(zhì)中受骨料的多重散射、反射的影響，造成的能量衰減。相對于均勻水泥漿介質(zhì)，混凝土中振幅不再線性衰減。

2.2 空洞模型

混凝土砌筑過程中，由于工藝技術(shù)、人為失誤及環(huán)境等因素的影響，常會有空洞缺陷的存在，為了研究混凝土中空洞對超聲傳播特性的影響，建立含空洞的混凝土模型。模型尺寸設(shè)置為400 mm×200 mm規(guī)格，骨料最大顆粒直徑為Amax=14 mm。空洞半徑為R0=20 mm(如圖7)。

圖6 抽取單道聲波波形記錄Fig.6 The single channel waveform record

圖7 含空洞混凝土模型Fig.7 Concrete model including round hole

圖8 波場快照對比圖Fig.8 The wavefield snapshot(a) t=28 μs時刻波場快照；(a) t=62 μs時刻波場快照

圖9 抽取單道聲波波形記錄Fig.9 The single channel waveform record

圖10 含裂隙混凝土模型Fig.10 Concrete model including horizontal cracks

圖8為超聲聲波在含空洞混泥土模型中傳播時不同時刻的波場快照。從圖8(a)可以看出，t=28 μs時刻，超聲縱波已經(jīng)完全傳播到空洞下空間，但因為混凝土與空洞兩者之間波阻抗相差很大，聲波在空洞界面發(fā)生全反射，在空洞周圍形成了明顯的反射波包絡(luò)面。同時根據(jù)惠更斯原理，超聲聲波在空洞界面上，以新的點源的形式繼續(xù)傳播，形成新的波前圓弧包絡(luò)面，但因為混凝土的衰減作用，這一部分的聲波能量降低，使得圓弧包絡(luò)面不連續(xù)，特別是空洞正下方的波形，由于全反射出現(xiàn)部分波形缺失。如圖8(b)所示：t=62 μs時刻，縱波初至波波前面?zhèn)鞑サ侥Ｐ拖陆缑?，被吸收層吸收，基本無反射波的形成；橫波初至波波前面?zhèn)鞑サ娇斩聪陆缑?，同樣，空洞正下方的波形由于受全反射的影響出現(xiàn)缺失；此時形成的橫波的全反射包絡(luò)面較明顯，而t=28 μs時刻形成的縱波的全反射波，由于能量的衰減，縱波的全反射包絡(luò)面不再明顯。為了更加清晰地分析波形傳播過程中能量的衰減情況，抽取了在含空洞混凝土模型橫向中心位置(200,50)處的單道波形(圖9)：超聲信號在t=22 μs時刻縱波初至波被換能器接收，波形振幅達(dá)到最大，在t=39 μs時刻在空洞上界面發(fā)生全反射，在t=47 μs時刻在空洞下界面形成反射波，但由于上界面的全反射作用，下界面反射波振幅較?。辉趖=54 μs時刻形成的波峰，是聲波橫波在空洞上界面形成的全反射形成的，此后，橫波在空洞中無法繼續(xù)向下傳播。相對于圖6混凝土模型中的單道信號：整體上超聲波能量也是隨著傳播距離逐漸衰減，但前期振幅減小較小，且振幅變化紊亂，是由于聲波在空洞上界面發(fā)生全反射，能量主要聚集在空洞上方，也造成波形傳播到空洞下方時，振幅整體表現(xiàn)較小?？斩创怪狈较蛏系姆秶笾驴梢酝ㄟ^兩波峰間的時間差確定。

圖11 波場快照對比圖Fig.11 The wavefield snapshot(a)t=26 μs時刻波場快照；(b)t=56 μs時刻波場快照

圖12 抽取單道聲波波形記錄Fig.12 The single channel waveform record

圖13 含傾斜裂隙混凝土模型Fig.13 Concrete model including tilted cracks

2.3 裂隙模型

在混凝土澆筑過程中，由于水泥砂漿凝固過程中的漲縮特性，混凝土在橫向和縱向上會存在應(yīng)力的作用，為了防止?jié)q縮力過大，通常會留有橫向和縱向的縫隙來緩沖漲縮力的作用。為了研究橫向裂隙和空洞對混凝土中聲波傳播影響的差異，模型尺寸設(shè)置為400 mm×200 mm規(guī)格，骨料最大顆粒直徑為Amax=14 mm。裂隙長度為50 mm，寬度為10 mm。裂隙的最小尺寸不能小于聲波最小波長(圖10)。

圖11是超聲聲波在含裂隙混凝土模型中傳播不同時刻的波場快照。從圖11(a)可與看出：超聲縱波傳播到裂隙位置處，由于裂隙與混凝土之間波阻抗差異較大，聲波會在裂隙上界面發(fā)生全反射，縱波形成的強反射波與橫波的初至波相互干擾，造成橫波波前的包絡(luò)面紊亂，同時由于裂隙上界面對縱波的強反射，造成裂隙正下方的縱波初至波包絡(luò)面缺失。圖11(b)顯示，超聲橫波傳播到裂隙后，在裂隙上界面同樣會發(fā)生反射，但橫波在裂隙中無法傳播，造成裂隙正下方的橫波包絡(luò)面波形缺失；同時因為裂隙橫向尺寸較大，在裂隙上界面發(fā)生全反射的范圍較廣，傳播到裂隙下方的聲波較少，形成一個能量較弱的梯形區(qū)域。為了更加清晰地分析波形傳播過程中能量的衰減情況，抽取了在含裂隙混凝土模型橫向中心位置(200,50)處的單道波形，如圖12所示：波形整體上能量隨傳播距離逐漸衰減，超聲信號在t=22 μs時刻的振幅波峰為縱波初至波，在t=32 μs時刻的振幅波峰為縱波在裂隙上界面全反射形成的反射波，由于裂隙寬度較小，縱波在裂隙中傳播的時間較短，在t=38 μs時刻形成的波峰即為縱波在裂隙下界面形成的反射波，隨后縱波振幅逐漸減??；在t=55 μs時刻出現(xiàn)的振幅較大，是橫波在裂隙上界面形成的強反射波。相對于圖9含空洞混凝土模型中抽取的單道信號，整體上振幅在前期衰減較慢，但在t=55 μs時刻之后，振幅迅速衰減，且減小值較大，這是由于縱波、橫波在裂隙界面形成全反射，且裂隙橫向截面較大，全反射范圍較大，使得裂隙上、下空間能量分布存在較大差異。

圖14 波場快照Fig.14 The wavefield snapshot(a)t=36 μs時刻波場快照；(b)t=72 μs時刻波場快照

圖15 抽取單道聲波波形記錄Fig.15 The single channel waveform record

2.4 傾斜裂隙模型

為了對比研究橫向裂隙和空洞對混凝土中聲波傳播影響的差異，設(shè)置傾斜裂隙模型，模型尺寸設(shè)置為400 mm×200 mm規(guī)格，骨料最大顆粒直徑為Amax=14 mm。傾斜裂隙長度為99 mm，寬度為9.66 mm。裂隙的最小尺寸不能小于聲波波長(圖13)。

圖14是超聲聲波在含傾斜裂隙混凝土模型中傳播不同時刻的波場快照。由圖14(a)可知：在t=36 μs時刻，超聲聲波傳播到傾斜裂隙位置處，超聲波在傾斜裂隙上界面發(fā)生強反射，反射波向背傾角的方向傳播，形成明顯的反射弧形包絡(luò)面，在傾斜裂隙下方，出現(xiàn)初至波包絡(luò)面的錯位、缺失，形成一處能量分布較小的區(qū)域，隨著波形的傳播，與裂隙的接觸面積越大，該區(qū)域面積越大。由圖14(b)可知：根據(jù)惠更斯原理，超聲在傾斜裂隙界面以新的點源的形式向外繼續(xù)傳播，t=36 μs時刻在傾斜裂隙下方形成的能量分布較小的區(qū)域，由于新的點源形成的波形的傳播，能量逐漸增大，逐漸形成比較完整、連續(xù)的包絡(luò)面。為了更加清晰地分析波形傳播過程中能量的衰減情況，抽取了在含傾斜裂隙混凝土模型橫向中心位置(200，50)處的單道波形，由圖15可知：超聲信號在t=22 μs時刻為縱波初至波，顯示的波形振幅最大；t=37 μs時刻為縱波在裂隙上界面發(fā)生強反射形成的反射波，t=44 μs時刻的波峰，推測是橫波初至波引起的。對比圖12中含橫向裂隙混凝土模型抽取的單道信號：傾斜裂隙由于傾角的存在，反射波會向背傾角的方向傳播，能量向背斜方向發(fā)散，在裂隙正上方(200,50)處接收到的能量較少，使得振幅衰減較大，振幅值較小。

3 結(jié)論

1)鑒于混凝土本身的多相非均勻特性，采用粘彈性方程進(jìn)行數(shù)值模擬。

2)基于彈性波方程，在粘彈性介質(zhì)中，添加粘性介質(zhì)中的拉梅系數(shù)，推導(dǎo)出粘彈性介質(zhì)中速度-應(yīng)力的有限差分表達(dá)形式。

3)通過數(shù)值模擬，研究聲波在混凝土內(nèi)傳播過程中的能量衰減情況，分析骨料配比對能量衰減的影響。

4)根據(jù)不同時刻的波場快照對比圖、能量分布及單道波形情況，分析空洞、裂隙、傾斜裂隙等不同缺陷對超聲傳播的影響，為推斷識別混凝土內(nèi)部不同結(jié)構(gòu)的缺陷提供依據(jù)。