侯婉婷， 符力耘*， 魏佳， 王志偉

1 中國石油大學(華東)深層油氣重點實驗室， 山東青島 266580 2 中國科學院油氣資源研究重點實驗室， 中國科學院地質與地球物理研究所， 北京 100029 3 中國科學院大學， 北京 100049

0 引言

高溫介質地震波傳播是地熱資源勘查、巖石圈高溫結構探測和深層-超深油氣勘探等領域的重要研究課題.一般而言，高溫引起巖石黏彈性，快縱波波速隨溫度增加而降低，其理論基礎是熱彈性理論，即溫度場彈性動力學.熱彈性動力學方程由位移場波動方程和溫度場擴散方程耦合構成.經典熱彈性波動方程(Biot，1956a；Deresiewicz，1957)基于傅里葉拋物型熱傳導定律，其熱平衡瞬間形成，忽略了溫度變化趨于平衡需要時間，導致熱彈介質中波的速度隨頻率無限變大的非物理現象(Savage，1966；Armstrong，1984).Lord和Shulman(1967)引入弛豫時間修正傅里葉熱傳導定律，簡稱L-S雙曲型熱輻射模型，解決了經典熱彈波動方程的非物理解問題.本文基于L-S雙曲型熱輻射模型，研究熱彈介質中彈性波傳播特征和各種熱參數對速度頻散和振幅衰減的影響.

溫度對巖石力學性質的影響以及熱彈介質中波傳播的研究由來已久.Ackerman等(1966)通過實驗驗證了Ward和Wilks(1951)提出的熱彈性介質中存在熱縱波(一種準靜態(tài)慢縱波).Rudgers(1990)對含弛豫時間的L-S型熱彈性動力學方程進行頻散分析，從理論上預測熱彈性介質中存在快縱波、橫波和熱波的傳播，兩個縱波為熱耗散衰減波而橫波不受介質熱特性的影響.Carcione等(2018)和Wang等(2019)分別利用偽譜法和格林函數法求解L-S雙曲型耦合熱彈波動方程，模擬熱彈介質中快縱波、熱波和橫波的傳播行為，指出了熱彈介質中的熱波與孔彈介質(Biot，1956b)中的Biot慢縱波具有相似的傳播特征.物理本質上，Biot波是由于孔隙介質中壓力驅動形成的流體局域流傳播行為；而熱波則是由于P波膨-縮傳播模式在介質中產生的溫度梯度，形成的局域熱擴散傳播行為.熱彈介質與孔彈介質的這兩種局域行為均引起彈性波能量的耗散和衰減.上述研究側重熱彈波動方程的物理解釋和數值求解方法，熱參數的影響分析比較簡單，僅給出了熱導率簡單變化的影響，證實模擬方法的有效性.

本文著重研究各種熱參數對熱彈性介質中彈性波傳播的影響，詳細分析了熱導率、熱膨脹系數、比熱、熱擴散系數等介質熱參數的大范圍變化(臨界變化)對彈性波頻散和衰減的影響.熱傳導是巖石圈最主要的熱量傳遞方式，其表征參數熱導率是巖石熱物理性質最重要的參數；熱膨脹系數表示巖石受熱膨脹的幅度，與巖石熱應力與熱應變密切相關；熱擴散率通過介質的比熱和熱導率計算確定.熱導率、熱膨脹系數和比熱構成了巖石的熱彈性基本系數.本文采用平面波頻散分析方法，研究了這些熱參數變化對熱彈介質中波速頻散與能量衰減的影響.然后利用熱彈性動力學頻率域的二階格林函數進行波場快照數值模擬，直觀展示熱參數對快縱波和熱波傳播特性的影響.本文研究為地熱資源勘查、巖石圈高溫結構探測和超深層油氣勘探提供理論依據，也為實驗室條件下兩類熱耗散衰減縱波的觀測提供參數依據.

1 熱彈性動力學方程

Biot(1956a)和Deresiewicz(1957)將彈性理論與熱傳導定律相結合，建立了經典的熱彈性理論，表征熱負荷作用下溫度變化與應力應變之間的耦合關系.對于各向同性熱彈性介質，其應變-位移關系和應力-應變關系定義如下：

(1)

和

(2)

式中，σij(i,j=x,y,z)、εij(i,j=x,y,z)和ui(i=x,y,z)分別為應力、應變和位移分量，T為相對于初始溫度T0的溫度增量，λ、μ為拉梅系數，fij(i,j=x,y,z)表示施加的外應力，γ=(3λ+2μ)αT，其中αT表示熱膨脹系數.可見與常規(guī)彈性介質本構關系不同，熱彈性介質中引入了溫度的影響，產生熱應力和熱應變，其大小取決于介質的熱膨脹系數.相應的動量守恒方程為：

(3)

式中，fi(i=x,y,z)表示體力分量，ρ表示體密度.

經典的傅里葉熱傳導定律可表示為：

(4)

式中，k表示熱導率，c表示恒應變下的比熱，q表示熱源，Δ為拉普拉斯算子.該定律精確描述了熱流與溫度差之間的關系，但由于熱流密度和溫度梯度同步，當介質產生溫度梯度時，無窮遠處瞬間感受到溫度變化，從而導致熱彈性介質中波以無限大速度傳播，即熱彈性波動方程存在非物理解，這在實際工程應用中問題比較突出(Qi and Suh，2010；Wu et al.，2013).

Lord和Shulman(1967)引入弛豫時間修正傅里葉熱傳導定律，簡稱L-S雙曲型熱輻射模型，解決了經典熱彈性波動方程的非物理解問題，對應的熱傳導方程可表示為：

(5)

式中τ表示弛豫時間，當τ=0時，(5)式退化為傅里葉熱傳導定律(4).將應力-應變關系帶入動量守恒方程，利用應變-位移關系，得到表征位移與溫度之間關系的熱彈性動力學方程：

(6)

而對于傳統(tǒng)拋物型傅里葉熱傳導定律，經典熱彈性動力學方程為：

(7)

由方程(5)可知，在熱傳導過程中，熱與應變耦合只涉及體應變，與剪切應變無關，建立的熱彈性動力學方程(6)僅考慮波與熱相互作用的一種機制：P波的膨-縮傳播模式在介質中產生溫度梯度，形成耗散衰減的慢縱波模式(T波)，這是實驗證實的一種熱與波作用機制.因此，在此理論框架下，剪切波(S波)不受溫度影響.

2 頻散分析

2.1 平面波頻散方程

為了分析熱彈介質中波的傳播特征，將位移矢量u與T的變化表示為平面波形式：

(8)

為簡單起見，若假定平面波為橫波，其位移矢量方向與傳播方向相互垂直，即Uili=0，熱彈方程的頻散關系簡化為(Deresiewicz，1957)：

μs2-ρω2=0，

(9)

式中復波數s=ω/vc，得到橫波的傳播速度：

(10)

由(10)式可知，橫波的傳播不受介質熱特性的影響.若假定平面波為縱波，其位移矢量方向與傳播方向一致，即Uili=1，得到縱波的頻散關系：

(11)

(11)式的解為：

(12)

(13)

對于經典熱彈性動力學方程，(12)式中M=iωa2，當ω→時，M→，由方程(10)可知會產生非物理解.而對于L-S廣義熱彈性理論，M=iωa2/(1+iωτ)，Rudgers(1990)在晶格模型中將弛豫時間表示為：

(14)

因此，當ω→時，M→a2/τ=，兩個縱波復波速為：

(15)

當角頻率ω=0時，由表達式(10)知，復速度為實數：

vc=vT=0，vc=vE=vA，

(16)

其中，vT表示熱波速度，vE表示快縱波速度.利用復速度vc，得到相速度vp與衰減因子A(Carcione, 2007)：

(17)

Deresiewicz(1957)引入衰減系數L，用來表示每個應力周期內，耗散的能量與總振動能量的比值：

(18)

2.2 平面波頻散分析

利用表1中的熱彈性介質模型參數進行頻散分析，通過2.1節(jié)平面波頻散方程可知，經典熱彈性波動方程以拋物型熱傳導方程為基礎，未考慮弛豫時間的影響，各向同性彈性連續(xù)體受到熱擾動時，在距熱擾動無窮遠處會立即感受到擾動，即速度無限大，其波速和衰減頻散曲線如圖1所示.可見，在熱效應作用頻段，快縱波隨頻率增加反而減小，而熱波速迅速增加至無窮大，其對應的衰減為非物理的倒弛豫峰.

將表1熱彈參數應用到雙曲型L-S熱彈理論得到的頻散曲線如圖2所示，可見，考慮加入熱擾動介質達到新穩(wěn)態(tài)的平衡時間(弛豫時間)后，熱效應作用頻段熱波速隨頻率增大趨于穩(wěn)定值，避免了經典熱彈性動力學中速度無限大問題.由圖2b的衰減頻散曲線可知，熱波在低頻時表現出強衰減的擴散傳播特性，而在高頻時表現為弱衰減的波動傳播特征.在高低頻之間熱效應過渡頻段及其頻散特征取決于熱介質的特征參數.

表1 熱彈性介質模型參數Table 1 Parameters for a thermoelastic medium model

圖1 經典熱彈波動方程的(a)速度頻散曲線和(b)振幅衰減曲線Fig.1 Velocity dispersion curves (a) and Amplitude attenuation curves (b) of classical thermoelastic wave equation

圖2 L-S雙曲型耦合熱彈波動方程的(a)速度頻散曲線和(b)振幅衰減曲線Fig.2 Velocity dispersion curves (a) and Amplitude attenuation curves (b) of L-S hyperbolic coupling thermoelastic equation

利用L-S方程分析熱彈介質中三種波的傳播，假設平面波入射，可知熱彈介質中的橫波不會受到介質溫度與熱彈參數的影響，只與介質密度及介質的彈性剪切模量有關，但由于實際介質非常復雜，橫波的傳播取決于很多因素，在本文所應用的條件下，橫波不發(fā)生頻散，當頻率約為100 MHz時，橫波與熱波傳播速度相同，如圖3所示.

圖3 L-S雙曲型耦合熱彈波動方程快縱波、熱波和橫波的速度頻散曲線Fig.3 The velocity dispersion curves of a fast longitudinal, thermal wave and shear waves of L-S hyperbolic coupling thermoelastic equation

2.3 熱彈參數影響分析

L-S熱彈理論中位移場與溫度場相互耦合，因此快縱波也受到溫度場的影響，快縱波和熱波的相速度和衰減系數根據平方根復速度(方程(17))計算得到，且由方程(12)可知其波速與衰減系數均涉及熱彈參數.基于表1中熱彈參數作為基準量，本文著重分析表征巖石熱物理性質的熱導率、熱膨脹系數及比熱對波傳播特征的影響.

熱導率表征巖石的導熱能力，不僅是研究地殼、上地幔熱結構及地球深部熱狀態(tài)的重要參數，在勘探開采等領域也得到廣泛應用(沈顯杰等，1988).熱導率是熱彈介質中波傳播最重要的影響參數，其大小不僅依賴于介質的礦物組成，還取決于礦物分布、幾何形狀、內部結構等因素.對于大多數巖石類型，熱導率隨溫度的升高而降低，與礦物結晶特性有關.對于孔隙介質和顆粒狀材料，熱導率還依賴于孔隙和晶粒尺寸等介質特征.通過經典熱彈理論與L-S熱彈理論對比可知，弛豫時間對熱彈介質中波的傳播性質有重要影響，通過分析與計算(Rudgers，1990)了解到熱導率與弛豫時間有很強的相關性.圖4描述了熱導率與快縱波和熱波的頻散與衰減的關系，數值計算過程中選取不同數量級的熱導率.可見，隨著熱導率的增加，相速度和衰減系數的頻散效應向低頻移動，但其低頻和高頻穩(wěn)態(tài)值并未發(fā)生改變，說明熱導率不改變波傳播的極限速度和衰減幅度；熱頻散過渡帶的頻寬也未見改變，但改變了熱頻散效應出現的頻率，介質熱導率越高，對應的熱頻散效應的初始頻率越低.需要強調的是熱波傳播特征發(fā)生了較大變化，從低頻慢速強衰減的擴散傳播模式跨越到高頻快速弱衰減的波動傳播模式，熱導率發(fā)揮了至關重要的作用.由方程(14)可知，弛豫時間與熱導率呈正相關，高熱導率情況下，介質受到熱擾動后達到穩(wěn)態(tài)的時間增大，即熱導率的升高增加了達到平衡的弛豫時間，從而降低了擴散頻率，如圖4b所示，衰減頻帶向低頻移動.因此，在實驗觀測中，高熱導率介質更容易觀測到熱波，下一節(jié)將通過格林函數數值模擬熱波傳播，直觀觀測熱導率的影響.

熱膨脹系數表征介質的熱膨脹量，不僅可以用來衡量介質的熱穩(wěn)定性能，還與介質受熱時所產生的熱應力與熱應變密切相關，因此熱膨脹系數對介質的結構有很大的影響，可能會導致結構破壞.與熱導率影響不同，熱膨脹系數對熱彈介質中波傳播特性的衰減影響主要表現在幅度方面.結合熱導率影響，可知，當介質熱導率較大時，可以觀測到波狀傳播的熱波.因此，取熱導率k=280 m·kg/(s3·K)對熱膨脹系數進行分析，如此大的熱導率值在一般巖石材料中是沒有的，在此僅作為參照值來闡述熱參數變化會對彈性波傳播特征產生多大的影響.圖5為不同熱膨脹系數快縱波和熱波的波速和衰減頻散曲線，與圖4熱導率影響特征相比，熱膨脹系數不改變相速度及衰減熱頻散效應發(fā)生的起始頻率及其過渡帶頻寬，說明熱膨脹系數并不影響擴散頻率，衰減頻帶未發(fā)生改變.由圖5a可知，熱膨脹系數改變了波傳播極限速度，特別是快縱波改變幅度較大.值得注意的是熱膨脹系數對兩縱波的影響方式相反，熱膨脹系數越大，快縱波相速度在全頻范圍內越大，而此時介質熱穩(wěn)定性能較低，熱波的相速度在高頻段越小.熱膨脹系數對熱波衰減系數影響很小，而對快縱波衰減系數的峰值影響稍大.

對于超高溫介質熱擴散率與比熱是不可或缺的參數，而熱擴散率可由熱導率和比熱計算得到.前面已經對熱導率和熱膨脹系數的波傳播影響進行了分析，本節(jié)研究比熱對熱彈介質中波傳播特性的影響，介質的比熱由介質組成成分決定.圖6為不同比熱介質快縱波和熱波的頻散曲線與衰減曲線，其中熱導率也為280 m·kg/(s3·K).可見，比熱的影響兼顧熱導率和熱膨脹系數的影響，既改變熱頻散效應發(fā)生的起始頻率，又影響極限速度的大小及快縱波弛豫峰值，但影響幅度要小得多.由圖6a可知，隨著比熱的增大，快縱波和熱波相速度的增量區(qū)均向高頻方向移動，但對兩縱波的影響相反.比熱越大快縱波的相速度在全頻范圍內越小，而熱波的相速度在高頻段越大；相比于快縱波，比熱對熱波相速度影響較小.且由圖可知，隨著比熱的增大，對縱波相速度的影響會越小.分析圖6b可知，隨著比熱增加，擴散頻帶向高頻移動，此時弛豫時間減小，達到穩(wěn)態(tài)的時間縮短.比熱綜合晶格振動和粒子熱運動的貢獻來表征巖石熱容性，是描述介質熱學性質的重要物理量，在常溫及高溫條件下，晶格振動影響占主導，因而對快縱波速度影響相對較大.

圖4 不同熱導率快縱波和熱波的(a)速度頻散曲線和(b)振幅衰減曲線Fig.4 Effect of thermal conductivity variety on Velocity dispersion curves (a) and Amplitude attenuation curves (b) of fast P and thermal waves

3 數值模擬

格林函數一般定義為偏微分方程自由空間中點源的解，是構成偏微分方程積分形式的基本解.針對L-S熱彈性動力學方程的格林函數研究鮮有發(fā)表(Wang et al.，2019)，其頻率域的二階格林函數張量解的簡要推導見附錄B，本文利用該格林函數張量進行熱介質波場快照數值模擬，直觀觀察熱彈性介質中縱波、橫波和熱波的傳播行為，傳播特征與第2節(jié)中的頻散分析結果一致.本文著重分析不同力源方向對三種波的波形特征及能量衰減的影響，同時再現熱彈介質參數對快縱波和熱波傳播特性的影響.

本節(jié)利用二階格林函數張量模擬波動傳播的快縱波和熱波，模型大小為NX×NZ=200×200，縱橫網格間隔均為0.1 mm，集中點源采用位于網格中點處主頻為3.5 MHz的雷克子波，時間步長dt=0.1 μs.

圖7表示在表1熱彈參數條件下，即熱導率為2.8 m·kg/(s3·K)，3 μs時刻的波場快照，(a)表示在水平力源條件下，水平方向位移的波場快照；(b)表示在垂直力源條件下，垂直方向位移的波場快照；(c)表示在熱源條件下，溫度場的波場快照.為了研究熱導率對快縱波和熱波傳播特征的影響，圖8中熱導率k=280 m·kg/(s3·K).由圖7和圖8可知，在力源作用下的波場快照中觀察到快縱波與橫波，在熱源作用下的溫度位移場中觀察到快縱波與熱波.由圖4可知，在該主頻下，當k=2.8 m·kg/(s3·K)時，熱波表現為慢速強衰減的擴散傳播模式，而k=280 m·kg/(s3·K)時，熱波以快速弱衰減的波動模式傳播，對于快縱波，高熱導率比低熱導率的傳播速度更快，而對于橫波，熱彈性介質參數不影響其傳播特性，均與格林函數數值模擬結果一致.

圖7 熱導率k=2.8 m·kg/(s3·K)，3 μs時刻的波場快照，震源分別為(a)水平力源，(b)垂直力源和(c)熱源Fig.7 Snapshots at 3 μs for k=2.8 m·kg/(s3·K), the source are a horizontal elastic force (a), a vertical elastic force (b) and a heat source (c), respectively

圖8 熱導率k=280 m·kg/(s3·K)，3 μs時刻的波場快照，震源分別為(a)水平力源，(b)垂直力源和(c)熱源Fig.8 Snapshots at 3 μs for k=280 m·kg/(s3·K), the source are a horizontal elastic force (a), a vertical elastic force (b) and a heat source (c), respectively

為了研究熱膨脹系數與比熱對快縱波和熱波傳播特性影響，在熱源條件下，熱導率k=280 m·kg/(s3·K)時，分別觀察介質熱膨脹系數不同和介質比熱不同時的波場快照.圖9a為熱膨脹系數αT=4.0×10-6K-1，對比圖8c和圖9a波場快照圖可知，熱膨脹系數對兩縱波的影響方式相反，在該主頻條件下，熱膨脹系數越大，快縱波速度越大，而熱波速度越小，與頻散分析結果(圖5)相符合.圖9b為比熱c=194 kg/(m·s2·K)，3 μs時刻的波場快照，通過與圖8c對比可知，比熱增大時快縱波速度減小，而熱波傳播速度增加，與圖6中的頻散分析結果一致.

圖9 熱源條件下(a)熱膨脹系數αT=4.0×10-6K-1和(b)比熱c=194 kg/(m·s2·K)時的波場快照Fig.9 Snapshots for c=194 kg/(m·s2·K) (a) and αT=4.0×10-6K-1 (b), where source is a heat source, respectively

4 討論與結論

地震波在熱彈介質中的傳播是巖石圈高溫結構探測、地熱資源勘查和深層-超深油氣勘探等領域的重要研究課題.基于傅里葉拋物型熱傳導定律的經典熱彈性動力學理論，由于忽略了溫度分布趨于平衡需要時間而導致波速隨頻率無限增大的非物理現象.基于L-S雙曲型熱傳導方程的熱彈性動力學理論，由于引入弛豫時間修正項而解決了這一非物理現象問題，本文通過對比這兩種熱彈性動力學理論的頻散曲線，分析了弛豫時間對熱彈介質中波傳播特征的影響.

含弛豫時間修正項的L-S雙曲型熱耦合熱彈波動方程，從理論上預測了在熱彈性介質中存在兩種熱耗散衰減縱波(即快縱波和熱波)和橫波，本文利用平面波頻散分析和數值模擬再現了這些波的傳播特征.熱波與孔彈介質中的Biot慢縱波(如，Ba et al.，2017，2019；Zhang et al.，2019)具有相似特征，兩種波分別為壓力和溫度的局部差異而產生的流體和熱流局域流，其頻散和衰減特征分別受制于介質的熱導率和流體黏度.兩種波低頻時呈現慢速高衰減的擴散特征，而高頻時表現為高速弱衰減的波動模式.本文的數值模擬表明，高熱導率材料在高頻實驗條件下可以觀測到熱波.Ackerman等(1966)在固體氦中觀測到熱波，McNelly等(1970)和Jackson等(1970)在NaF晶體中觀察到熱波.作為熱彈性問題不可分割的理論組成部分，熱波的客觀存在具有重要的科學意義(Gueguen，2013).

鑒于熱彈參數對熱彈介質中彈性波傳播的重要影響，本文詳細分析了熱導率、熱膨脹系數和比熱三個基本參數的變化對波速頻散和振幅衰減的影響.結果表明，熱導率決定了兩種熱耗散縱波波速頻散與能量衰減的臨界變化；熱膨脹系數主要影響快縱波的波速和熱波衰減的幅度，而比熱既影響波速和衰減的臨界變化，又影響變化的幅度.格林函數方法數值模擬展現了熱彈介質中橫波、快縱波和熱波的波場快照，熱參數(熱導率、熱膨脹系數和比熱)對快縱波和熱波傳播特征的影響與平面波頻散分析結果一致.本文研究為地熱資源勘查、巖石圈高溫結構探測和超深層油氣勘探提供理論依據，也為實驗室條件下兩類熱耗散衰減縱波的觀測提供參數依據.

附錄A 文中涉及的符號列表

符號含義符號含義ui位移分量(i=1,2,3)ω角頻率εij應變分量(i, j=1,2,3)vc復速度σij應力分量(i, j=1,2,3)d慢度矢量λ, μ拉梅系數x位置矢量

續(xù)表

附錄B

關于熱彈問題的格林函數，前人做了大量研究，主要涉及熱靜力學問題格林函數(Kupradze et al.，1976)和時間域及頻率域的經典熱彈性耦合熱彈波動方程的格林函數(Nowacki,1975；Tosaka,1985；Tosaka and Suh,1991；Norris,1994；Scarpetta et al.，2014)，未考慮弛豫時間，存在非物理解問題. Wang等(2019)推導了含弛豫時間修正項的L-S雙曲型熱耦合熱彈波動方程，得到了均勻各向同性材料中點載荷(力或熱源)對應的二階張量格林函數.

若不考慮熱流的影響，方程(6)簡化為：

(B1)

對耦合系統(tǒng)運用傅里葉變換進行求解，傅氏變換為：

(B2)

得到熱彈性動力學方程在頻率域的雙曲型微分控制方程(為了方便，在下面的公式表達中去掉-)：

(B3)

將(B3)式寫成矩陣方程形式：LijUj=0，考慮二維情況，經過數學計算得到系數算子Lij：

其中Di=?/?xi(i=1,2).

(B4)

式中l(wèi)ij為Lij的代數余子式，因此系數lij為：

(B5)

式中βij為lij的余子式.將(B5)式代入方程(20)得到：

ΛΦ*=-δ(x-y)，

(B6)

式中微分算子Λ為：

(B7)

(B8)

其中：

(B9)

(B10)

其中：

(B11)

K0，K1，K2分別為修正的零階、一階和二階貝塞爾函數，并且：

(B12)