利用有限元熱應力計算方法對熱事件月震成因的探討*

張君策,胡才博,石耀霖

(中國科學院大學地球與行星科學學院 中國科學院計算地球動力學重點實驗室, 北京 100049) (2021年10月14日收稿; 2022年4月22日收修改稿)

月球作為地球的衛星,一直受到各界科學家的關注。由于航天技術的進步,人們已經實現多次登月活動。在美國Apollo多次探月計劃中,人類通過各種技術手段獲取月球上空及表面的觀測數據,采集月壤和月巖樣品,有助于人們研究月球表面及其內部溫度的時空演化特征,進而研究與月球熱狀態有關的月震活動性。

阿波羅登月計劃對月球的重力場、電場、磁場進行了全面探測,并布設相當數量的月震儀器得到了相對全面的月震觀測數據。根據前人對Apollo12、14、15、16、17的研究統計,截至2018年,長周期月震儀共探測到13 058次月震事件,其中深源月震7 083次(深度大約在700～1 200 km范圍),深度大概幾十到200 km范圍內的月震28次,隕石撞擊1 743次,熱事件555次[1-2]。Duennebier和Sutton[3]認為,月表熱事件多數震級較小、震源較淺,這些月震是由于月表日周期溫度變化導致的熱應力引起的熱月震。有學者利用新的模式識別算法對Apollo17號記錄的時長為8.3個月的月震事件進行分析,發現它們的發生與日出日落時間有一定的相關性,認為上述月震事件的物理解釋可能是熱驅動[3-4]。

一般認為,由于月球淺表受到太陽輻射周期性變化導致溫度升降和熱應力擠壓-拉伸的周期性變化,從而導致月球發生熱事件月震[3,5]。因此,熱事件月震的研究需要對月球表面溫度的時空分布有更定量的了解。Apollo15的鉆孔溫度探測數據表明,月表該處最高溫度為374 K,最低溫度為92 K,而Apollo17著陸點測得月表最高溫度為384 K,最低溫度為102 K[6-7]。由于太陽輻射作用,月球表面溫度與緯度相關,赤道的溫度在滿月時大概為390 K,而月夜期間,月球不受到太陽輻射,溫度下降到100 K左右[8-9]。

月球表面大氣稀薄,月表溫度場基本上不受空氣對流的影響?？傮w而言,月球的溫度場主要有3個能量來源：一個是白天太陽光照的熱量,一個是來自月幔的熱量,最后一個就是地球對月球輻射所產生的熱量輸入[10]。對于月球表面而言,溫度場的時空變化基本不受來自月幔的熱量影響,而起重要影響的因素就是向內吸收的太陽輻射和向外釋放的黑體輻射。月球的自轉存在29.5 d的自轉周期,在白晝條件下,月球表面的總熱量來源主要為太陽光照與月表黑體輻射的熱量之差,在夜晚條件下,月表只存在月表向外的黑體輻射。月球表面接收的太陽輻射量呈周期性變化,因此,月表溫度也隨著月球日周期性變化。

月球不同區域表面存在不同厚度的月壤,而月壤的熱導率遠低于月巖,月壤厚度對月球表面附近的溫度的時空分布影響很大。Mitchell和Pater[11]分析了不同深度下的溫度變化因素,溫度高于350 K時,熱輻射對月球淺層溫度變化影響起主導作用,而在下層的固體巖石中,熱傳導起主導作用。在此基礎上,Vasavada等[12]建立了對應的雙層模型,認為月球表面溫度在白晝時主要依賴太陽光照,日落后,月球表面失去熱量來源,月表溫度迅速下降。綜合目前的研究,月球上大部分區域的月壤厚度約為幾米,在月海區一般為3～5 m,而月球高地月壤厚度可能達到10 m[13-14]。Hayne等[15]對月壤特征厚度分布進行統計分析,得出月壤特征厚度的分布呈現非對稱性的概率分布,在0～3 cm特征厚度范圍內的月壤區域相對較少,而月壤特征厚度的最概然值約為7 cm。因此可以認為,月壤平均特征厚度可選取5 cm作為典型值。

目前針對月表溫度變化的研究多為基于均勻介質模型進行分析,對溫度場的演化進行研究[12,16-17],或是基于Apollo樣本以及近似模擬材料,對月球表面月壤風化層性質進行探究[15,18]。而對溫度變化引起的熱應力變化更是缺乏深入的定量研究。實際上,月壤材料相對于月球基巖,熱力學性質存在相當大的差異。根據月球波速結構的研究結果,月壤相對于月球基巖來說,月表介質的彈性波波速僅為深部的1/10[19],而月壤相對松散,孔隙度較高,因而導熱率只有基巖導熱率的1/1 000[10]。月表不同區域月壤厚度存在明顯的差異[15,17-18]。因此,研究月表溫度場以及熱驅動的熱月震,需要對月壤的特征厚度和熱力學參數有更定量的認識。本文利用自主研發的二維熱彈性耦合有限元程序,以溫度變化幅度最大的月球赤道區域為例,定量研究月球淺表的溫度、垂直位移及熱應力的時空分布,分析月壤特征厚度對溫度、位移及熱應力演化的影響,進而探討月表附近熱事件的動力學成因。

1 模型與方法

本文不考慮來自月幔的熱量,只考慮月表向內吸收的太陽光照熱量和向外釋放的黑體輻射熱量。根據觀測和模擬計算,月表溫度變化會呈現出29.5 d的周期性[16,20]。建立二維熱彈性耦合有限元模型,考慮月壤的不同特征厚度及熱力學參數隨深度、溫度變化的非線性,研究4個完整的月球日月表月壤的溫度、位移和熱應力的時空變化。相關的熱力學參數的理論公式如下。

Hayne等[15]根據行星數據系統的研究結果,給出月壤密度模型

(1)

其中：ρ(z)表示隨深度變化的月壤密度,ρd表示在月壤底部處的密度,ρs表示月壤表面密度,z表示深度,H表示月壤特征厚度。

Schreiner等[18]利用Apollo計劃帶回的樣本數據,對月表風化層的性質做了回歸分析,并給出比熱的擬合公式

Cp(T)=c0+c1T+c2T2+…+cNTN.

(2)

其中：Cp(T)表示熱容隨溫度變化函數,T表示溫度,c0,c1,…,cN表示擬合公式的比熱系數。

根據實驗數據,熱導率與溫度密度存在相關函數關系,得出以下公式[11,15]：

(3)

(4)

其中：κ為月壤的熱導率,κc是與深度相關的熱導率系數,κs和κd分別表示月壤表面和底部的熱導率,ρs和ρd表示月壤表面和底部的密度,χ表示輻射熱導率系數,T為月壤溫度。

本文中,ρs=1 100 kg/m3,ρd=1 800 kg/m3,κs=7.4×10-4W·m-1·K-1,κd=3.4×10-3W·m-1·K-1,χ=2.7[17]。

月表太陽吸收熱量的計算公式[21]如下

QS=(1-A)F(t).

(5)

其中：QS為月表吸收太陽光照的熱量,A為月表太陽光反射率,F(t)為月表太陽光入射熱量。

本文選取月表太陽反照率經驗公式[15,22-23],此公式已得到檢驗：

(6)

其中：A(θ)表示月表太陽反射率,θ表示太陽入射角,以弧度為單位;A0為反照率擬合基數;a、b均為擬合系數。本文采用A0=0.12,a=0.06,b=0.25[17]。

月表太陽光入射熱量F(t)的計算公式[17]如下

F(t)=

(7)

其中：S0為太陽常數,RAU表示太陽和月球間的距離,Φ表示緯度,δ為太陽偏角,h為與時間有關的時間角,h=2πt/P,t為時間,P為月球日周期。本文中以赤道地區為例,Φ=0°,δ=0°,S0=1 361 W/m2,RAU=1.0[17]。

本文使用Stefan-Boltzmann定律描述月球表面向外輻射的熱量QT

QT=εσT4.

(8)

其中：ε表示物體發射率,σ表示Stefan-Boltzmann常數,T表示溫度。本文中,ε=0.95,σ=5.67×10-8W/(m2·K4)[17]。

熱應力是由于物體在存在溫度差的時候,溫度的差異性分布會導致物體內部存在一定的局部性膨脹或收縮現象,溫度的平衡過程會產生一定的熱應力。對于月壤彈性體而言,溫度演化公式[17]如下

(9)

其中：ρ,Cp,κ,Q分別為月壤的密度、比熱、熱導率和熱生成率。

月壤彈性體的應力平衡方程為

σij,j+fi=0.

(10)

彈性本構關系為

σ=[D](ε-εα).

(11)

其中：σ和ε分別表示應力和應變的列向量形式,εα表示熱應變,[D]為彈性矩陣。彈性矩陣表示方法如下

(12)

其中：E,ν,G分別為楊氏模量、泊松比和剪切模量。

熱應變可以表示為

(13)

其中：α為熱膨脹系數,ΔT為相對于初始溫度的溫度變化。

所以,二維平面應變熱應力本構方程如下

(14)

在有限元計算中,邊界條件是計算的一個重要因素。本文存在3個物理場,其一為溫度場,左邊界、右邊界、下邊界為絕熱邊界,上邊界給定邊界的熱流邊界條件,具體表示如下

Q0=QS-QT.

(15)

其中：QS、QT的表達式見式(5)和式(8)。

另2個物理場為位移場、應力場,其中左右邊界水平方向上位移約束為0,垂直方向上不做約束;下邊界垂直方向上位移約束為0,水平方向位移不做約束;上表面自由。

2 有限元模型的設計

為探究月表附近熱事件的動力學機制,本文采用熱彈性有限元程序,建立了二維月壤剖面模型,水平0.1 m,深度1 m,水平網格尺寸為2 mm,0.2 m以淺的垂直網格尺寸為1 mm,0.2 m以下的垂直網格尺寸為4 mm。采用結構化四邊形網格,共有20 451個有限元節點,20 000個四邊形單元。分別設計了4個不同的有限元模型,其月壤特征厚度H分別為0、2、5、7 cm,其中月壤特征厚度為零的模型表示密度ρ取為底部密度ρd、熱傳導率κc取為底部熱導率κd。其余的幾何、材料、初始條件、邊界條件一致。初始溫度的分布情況對月球淺表溫度場數值模擬的影響很大。廖翼[24]研究不同的初始溫度分布(均勻條件或者指數條件)對月表附近溫度場模擬計算的收斂性的影響,發現采用300 K的均勻初始溫度需要120個月球自轉周期才能達到計算穩定。本文采用月壤均勻初始溫度為255 K,模型底部溫度為255 K,表面采用月表太陽吸收熱量Qs和月球表面向外輻射的熱量QT相疊加的熱流邊界條件,兩個側面是無熱流邊界條件。對于位移場,月表自由,兩個側面和底面法向位移約束為零,切向位移自由。應力的計算采用式(11)。

本文材料模型選取見表1。

表1 月壤材料參數Table 1 Lunar soil material parameters

一個完整的月球日約為29.5個地球日。為了保證非線性計算收斂性以及耦合場計算精度,本文選取177 s為計算步長,相當于月球的0.1 min,每60個計算時間步存儲1次計算結果,相當于月球上的每6 min存儲1次計算結果。4個有限元模型的起算時間均為月球上的正午12：00,總共的計算時間為1.019 52×107s,相當于計算4個完整的月球日。為提高計算效率,采用并行程序,利用9個計算核數,計算不同月壤特征厚度的4個模型的溫度、位移和熱應力場。

3 模擬結果

選取月壤特征厚度為2 cm的有限元模型為參考模型,首先給出參考模型的溫度、位移和熱應力的計算結果。在太陽熱輻射吸熱和月球本身熱輻射放熱的共同作用下,月球表面及淺表的溫度、位移和應力都呈現周期性的變化。對于我們的二維模型而言,深度方向是1 m,在垂直方向上,月表是自由的,月壤因溫度變化將自由伸縮,沒有約束力,因此理論上,垂直正應力為零。由于月表非線性的熱流邊界條件,有限元模擬得到的垂直正應力在月表附近可能會存在一定的數值誤差,并隨深度迅速衰減到零,故本文不考慮垂直正應力的影響。

圖1為深度分別為0,0.1,0.5,1,5,10,50 cm溫度在4個完整的月球日和第4個月球日內隨時間的變化曲線圖。圖1(a)和1(b)的起始時間是正午12：00(下同)。除第1個月球日受初始場影響以外,其余3個月球日,月表及不同深度的溫度,均呈現一定的周期性變化(圖1(a))。月表(圖1中的紅線)在1個完整的月球日(以第4個完整的月球日為例,橫坐標為72～96月球小時,其中橫坐標72和96剛好對應月球上的正午12：00),正午12點溫度達最大值384.85 K,早上6點溫度達最小值100.13 K。從正午12：00到第2天的早上06：00,月表溫度隨時間都是降低的,其中白天的12：00(橫坐標72)到18：00(橫坐標78),溫度從極大值開始急劇下降,18：00(橫坐標為78)到第2天早上06：00(橫坐標為90),溫度緩慢下降,到06：00達到極小值。從第2天早上06：00(橫坐標為90)到第2天正午12：00(橫坐標為96),溫度從極小值急劇升高到極大值,此過程與之前的12：00到18：00的過程剛好相反(圖1(b))。

圖1 2 cm特征厚度月壤參考模型,不同深度下月壤溫度隨時間變化曲線圖Fig.1 Curves of lunar soil temperature with time at different depths (Reference model of lunar soil with characteristic thickness of 2 cm)

月表之下0.1,0.5,1,5,10,50 cm的溫度隨時間的演化曲線,與月表的有一定的相似性,深度越大,溫度曲線的變化幅度逐漸指數式減少。深度為50 cm時,溫度基本不隨時間變化。與此同時,溫度的極大值出現比地表位相滯后,深度越深,極大值越往后推遲,0.1,0.5,1,5,10 cm的溫度極大值和出現時間分別為：(0.1 cm, 380.53 K,12：12)、(0.5 cm,367.75 K,12：30)、(1 cm,356.63 K,12：54)、(5 cm,310.33 K,15：00)、(10 cm,281.09 K,17：30)。溫度的極小值在1 cm以前依然是在06：00出現,5 cm以下極小值不在06：00出現,而是有所延遲(5 cm,195.70 K,06：36)、(10 cm,226.21 K,08：18)(圖1(b))。

Turcotte和Schubert[29]給出了地表單一頻率周期性隨時間余弦分布的溫度邊界引起的溫度幅值隨深度的指數衰減和相位隨深度的延遲的解析解,表明溫度幅值隨深度指數式減小,周期大的低頻成分衰減較慢,周期小的高頻成分衰減較快。本文的月表溫度含有多種頻率成分,因此,其溫度隨深度的衰減和相位延遲比單一頻率的解析解結果要更復雜一些。

圖2給出不同深度,參考模型的垂直位移隨時間或溫度的變化曲線,同樣呈現周期性變化(圖2(a))。垂直位移大于零代表抬升,小于零代表下降。以模擬的第4個月球日(同上)為例,月表的垂直位移變化范圍最大,深度越深,垂直位移的變化范圍越小,50 cm深度處的垂直位移幾乎為零(圖2(b))。月表垂直位移極大值也呈現位相差異(圖2(b))。經過1個月球日的熱脹冷縮,不同深度的垂直位移和溫度的曲線基本是大小不一的封閉閉環(圖2(c)),由于月表總位移是不同深度不同溫度熱脹冷縮變化量的積分求和,所以呈現的深度衰減和位相滯后形式更為復雜。

圖2 2 cm特征厚度月壤參考模型,不同深度下月壤垂直位移隨時間或溫度的曲線圖Fig.2 Curves of vertical displacement of lunar soil with time or temperature at different depths (Reference model of lunar soil with characteristic thickness of 2 cm)

圖3給出不同深度,參考模型的水平正應力隨時間或溫度的變化曲線,同樣呈現周期性變化(圖3(a))。水平正應力大于零表示拉張,小于零代表擠壓。以模擬的第4個地球日(同上)為例(圖3(b)),月表(紅線)的水平正應力的變化范圍最大,極值及時刻如下：(0 cm,2.066 MPa,06：00),(0 cm,-1.731 MPa,12：00),(0.1 cm,1.982 MPa,06：00),(0.1 cm,-1.674 MPa,12：12),(0.5 cm,1.742 MPa,06：00),(0.5 cm,-1.503 MPa,12：30),(1 cm,1.541 MPa,06：00),(1 cm,-1.355 MPa,12：54),(5 cm,0.791 MPa,06：36),(5 cm,-0.738 MPa,15：00),(10 cm,0.383 MPa,08：18),(10 cm,-0.348 MPa,17：30),其中正值表示拉應力極大值,負值為壓應力的極值(圖3(b))。月表水平正應力拉升的最大值出現在06：00,值為2.066 MPa;月表水平正應力擠壓的最大值出現在12：00,值為-1.731 MPa。0～0.5 m深度,水平正應力也是周期性拉壓的(圖3(b))。值得注意的是,18：00左右,水平正應力的變化最大,由壓應力迅速地變為較大的壓應力,然后經過晚上的12 h,拉應力緩慢增長到早上06：00的極大值(圖3(b))。經過1個月球日的熱脹冷縮,不同深度的水平正應力和溫度的曲線基本是長短不一的線段,月表對應的線段最長,50 cm深度對應的線段幾乎縮為1個點。高溫代表壓縮,低溫代表拉張,兩者是線性關系,溫度可以作為水平正應力的指示器(圖3(c))。如果考慮月表附近的拉應力破壞的話,18：00左右是拉應力增長最快的時刻,而早上06：00左右是拉應力最大的時刻,都是值得關注的。

圖3 2 cm特征厚度月壤參考模型,不同深度下月壤水平正應力σxx在4個月球日隨時間變化曲線圖Fig.3 Curves of horizontal normal stress σxx of lunar soil with time in four lunar days at different depths(Reference model of lunar soil with characteristic thickness of 2 cm)

圖4為2 cm厚月壤參考模型模擬的第4個月球日(同上)不同時刻(每2個月球小時)溫度隨深度變化曲線圖?？梢钥闯鰷囟戎饕?0 cm以淺的深度劇烈變化,變化范圍約為100.13～384.85 K,50 cm深度以下溫度基本不受月表影響,保持不變,約為255 K。綜合來看,08：00到16：00是白天,月表有太陽光照射接收熱量,溫度較高,其中12：00的月表溫度最高。18：00到06：00是夜晚,月表沒有太陽光照射不接收熱量,溫度較低,其中18：00左右溫度下降最快,06：00的月表溫度最低。

圖4 2 cm特征厚度月壤參考模型,不同時間溫度隨垂直坐標的變化曲線圖Fig.4 Curves of temperature with vertical coordinates for different time periods (Reference model of 2 cm characteristic thickness lunar soil)

圖5為2 cm厚月壤參考模型模擬的第4個月球日(同上)不同時刻(每2個月球小時)水平正應力隨深度變化曲線圖?？梢钥闯鏊秸龖χ饕?0 cm以淺的深度劇烈變化,水平正應力的變化范圍約為-1.731～2.066 MPa,既有拉應力,也有壓應力,拉應力的極大值要大于壓應力的極大值。50 cm深度以下水平正應力基本不受月表的影響,幾乎為零。綜合來看,08：00到16：00是白天,月表有太陽光照射接收熱量,月表水平正應力是壓應力,其中12：00的月表的水平壓應力最大(圖5(a))。18：00到次日06：00是夜晚,月表沒有太陽光照射不接收熱量,月表水平正應力是拉應力,其中18：00左右拉應力增長最快,06：00的拉應力最大(圖5(b))。若將垂向坐標軸改為相應的溫度軸,則無論是白天,還是夜晚,水平正應力與溫度都是在一條直線上,兩者是線性關系,溫度可以作為水平正應力的指示器(圖5(c)、5(d))。

圖5 2 cm特征厚度月壤參考模型,不同時間水平正應力σxx隨垂直坐標或溫度的變化曲線圖Fig.5 Curves of horizontal normal stress σxx with vertical coordinate or temperature at different time(Reference model of lunar soil with characteristic thickness of 2 cm)

圖6為2 cm厚月壤參考模型模擬的第4個月球日不同時刻(每2個月球小時)垂直位移隨深度變化曲線圖?？梢钥闯龃怪蔽灰浦饕?0 cm以淺的深度劇烈變化,變化范圍約為-0.017 0～0.012 0 m,既有抬升,也有下降,下降的極大值要大于抬升的極大值。50 cm深度以下垂直位移基本不受月表的影響,幾乎為零。注意,與溫度、水平正應力不同,月表的垂直位移的抬升、下降并不是跟白天和夜晚嚴格對應,而是有一定的延遲或滯后。從正午12：00,抬升量呈(12：00,0.006 3 m),(14：00,0.011 0 m),(16：00,0.011 6 m),(18：00,0.005 6 m)變化,均表示地表抬升(圖6(a))。而(20：00, -5.236 6e-04 m), (22：00, -0.004 9 m),(00：00,-0.008 5 m),(02：00,-0.011 7 m),(04：00,-0.014 6 m),(06：00,-0.017 3 m),(08：00,-0.010 3 m),(10：00,-0.001 4 m)變化,均表示地表下降(圖6(a)、6(b))。若將垂向坐標軸改為相應的溫度軸,則無論是白天,還是夜晚,垂直位移與溫度呈發散狀分布,兩者不是線性關系(圖6(c)、6(d))。

圖6 2 cm特征厚度月壤參考模型,不同時間垂直位移v隨深度變化曲線圖Fig.6 Curves of vertical displacement v with depth at different time (Reference model of lunar soil with characteristic thickness of 2 cm)

圖7給出不同月壤特征厚度對月表溫度演化的影響?？梢钥闯?白天12：00—18：00(圖7橫坐標72～78),06：00—12：00(圖7橫坐標90～96),僅在06：00和18：00附近有較大差異,除此之外,在白天的時段,4條溫度曲線幾乎完全重合。而夜晚時段18：00—06：00(圖7橫坐標78～90),月球接收不到太陽光照,只有月球向外輻射熱量,而此時,不同厚度的風化層保溫效應不同,也因此表現為不同的月壤特征厚度對溫度的時間演化有影響。以模擬的第4個月球日(同上)為例,夜晚時段,在早上06：00時,4個模型的溫度差異最大,分別為(0 cm,103.93 K),(2 cm,100.13 K),(5 cm,95.55 K),(7 cm,93.83 K)。

圖7 第4個月球日不同月壤特征厚度模型的月球表面溫度隨時間變化曲線(起始時間是12：00)Fig.7 On the fourth lunar day, curves of lunar surface temperature with time under different lunar soil characteristic thickness models (the starting time is 12：00 at noon)

圖8給出不同月壤特征厚度對月表水平正應力演化的影響?？梢钥闯?白天僅在06：00和18：00附近水平正應力有較大差異,除此之外,在白天的時段,4條曲線幾乎完全重合。而夜晚時段18：00—06：00(圖7橫坐標78～90),水平正應力與溫度類似,與月壤風化層厚度密切相關(圖8(a)、圖8(b))。

圖8 不同月壤特征厚度模型的月球表面水平正應力隨時間變化曲線(起始時間是12：00)Fig.8 Time varying curves of horizontal normal stress on the lunar surface of different lunar soil characteristic thickness models (the starting time is 12：00 at noon)

圖9、10、11分別給出不同月壤特征厚度模型在月球時間12：00、18：00、06：00時溫度、水平正應力和垂直位移及其與零特征厚度模型差異隨深度變化曲線?？梢钥闯?不同特征厚度模型的溫度、水平正應力和垂直位移隨深度都有較為明顯的變化,尤其是在較淺的月表附近。

圖9 不同月壤特征厚度模型的溫度隨深度的變化Fig.9 Curves of temperature with depth of different lunar soil characteristic thickness models

在正午12：00,10 cm以淺的區域,月壤特征厚度越大,相對零特征厚度模型,溫度越低,在5 cm深度附近差異可達-30 K(圖9(a)、圖9(b)),水平正應力的差異越大,在5 cm深度附近差異可達+0.4 MPa(圖10(a)、圖10(b)),而在10 cm以下的深度,溫度、水平正應力呈現相反的變化趨勢。在正午12：00,月壤特征厚度越大,相對零特征厚度模型,垂直位移增長越大,在10 cm深度處的最大差異接近+10 mm(圖11(a)、圖11(b))。

圖10 不同月壤特征厚度模型的水平正應力隨深度的變化Fig.10 Curves of horizontal normal stress with depth of different lunar soil characteristic thickness models

圖11 不同月壤特征厚度模型的垂直位移隨深度的變化Fig.11 Variation curve of vertical displacement with depth by different models

在傍晚18：00,10 cm以淺的區域,月壤特征厚度越大,相對零特征厚度模型,溫度越低,在5 cm深度附近溫度差異量可達+20 K(圖9(c)、圖9(d)),水平正應力的差異越大,在5 cm深度附近差異可達-0.3 MPa左右(圖10(c)、圖10(d)),而在10 cm以下的深度,溫度、水平正應力呈現相反的變化趨勢。在傍晚18：00,月壤特征厚度越大,相對零特征厚度模型,垂直位移增長越大,在20 cm深度處的最大差異接近+6 mm(圖11(c)、圖11(d))。

在早上06：00,10 cm以淺的區域,月壤特征厚度越大,相對零特征厚度模型,溫度越低,在5 cm深度附近溫度差異量可達+30 K(圖9(e)、圖9(f)),水平正應力的差異越大,在5 cm深度附近差異可達-0.4 MPa(圖10(e)、圖10(f))。在早上06：00,月壤特征厚度越大,相對零特征厚度模型,垂直位移增長越大,在月表處的最大差異接近+14 mm(圖11(e)、圖11(f))。

4 討論

1)2 cm月壤特征厚度的參考模型計算結果,月表溫度最大值為384.85 K,剛好是12：00,溫度最小值為100.14 K,剛好是06：00,代表的是月表赤道附近的溫度范圍。Apollo15的鉆孔溫度探測數據表明,月表該處最高溫度為374 K,最低溫度為92 K,而Apollo17著陸點測得月表最高溫度為384 K,最低溫度為102 K[6-7]。赤道的溫度在滿月時溫度大概為390 K,而月夜期間,月球不受到太陽輻射,溫度下降到100 K左右[8-9]。我們數值模擬的結果與觀測結果基本一致。

2)月壤上下邊界存在一定的溫度差,且由碎屑為構成成分的月壤孔隙度高且變形能力相對下層基巖而言較弱,因而使得月壤層內本身應力差較小。而如果月壤厚度較大,月表熱量難以及時傳導到月球深部基巖,進而使得巖石溫度變化較小。對應在熱應力中,則表現為難以產生較大的熱應力值,所以月壤較厚的區域熱事件分布相對稀疏,數量較少;而對于較薄的月壤,溫度相對而言易于傳播,使得月球基巖的溫度變化劇烈,產生較大的熱應力,因而薄月壤區域容易發生較為密集的熱事件[3]。月表溫度急劇下降最快的時刻是18：00左右,對應的水平正應力是拉應力,而且增加最快。月表溫度最低的時刻是06：00左右,對應的水平正應力是拉應力,而且此時的拉應力最大。18：00到第2天的06：00,水平正應力都是拉應力,而且是緩慢增長到最大的拉應力。因此,需要關注早晚6點左右的時刻,比較有利于發生月表附近的熱事件。

3)太陽照射的時段06：00—18：00,而且吸收熱量與太陽的入射角有關,12：00的入射角為0°,此時月表溫度達到極大值,06：00—12：00,入射角從90°逐漸變為0°,吸收的熱量越來越大,溫度迅速升高;12：00—18：00,入射角從0°逐漸變為90°,吸收的熱量越來越小,溫度迅速降低。18：00—06：00,夜晚,沒有太陽照射,完全是月球向外的熱輻射,溫度繼續降低,但是降溫的速度比較緩慢。隨著深度增加,溫度差異逐漸減小,最終達到平均溫度。對于不同厚度的月壤模型,月壤厚度越大對應的月表和月巖溫度差異越大,這與前人的研究結果相符[25,30]。這也與月球月壤的主要構成成分為碎屑相關,相對月巖來說導熱能力較弱。

4)月表垂直位移與月表溫度、水平正應力的時間曲線有所不同,有一定的滯后性。換句話說,月表垂直位移相對于溫度、水平正應力而言,有一個積分響應,有一個累積的過程。不同的月壤特征厚度,對于溫度和正應力而言,在白天影響較小,夜晚有一定影響,06：00左右4個模型的差異最大。不同的月壤特征厚度,對于垂直位移而言,4個模型的差異對所有時段都有顯著差別,更是體現垂直位移的時間滯后性和累積性。

5)本文在月球淺表溫度計算中沒有考慮來自月幔的熱流,因為50 cm以下深度的月熱流密度實測值僅為20～40 mW/m2[31],而計算(第4個月球日)表明,08：00從月表向下傳導的1 m深度平均熱流密度為39 mW/m2,18：00從底部向月表傳導的1 m深度平均熱流密度可達87 mW/m2,月表1 m深度范圍內的平均熱流差異也很大,甚至反向,模擬得到的熱流密度與實測結果在數量級上一致。模擬結果還表明,08：00從月表向下傳導的1 cm深度平均熱流密度達3.7 W/m2,18：00從1 cm深度向月表傳導的平均熱流密度達7.2 W/m2,遠遠超過50 cm以下深度的月熱流密度實測值20～40 mW/m2,說明來自月幔的熱流對于月表附近的溫度演化可以忽略不計[25]。

6)本文建立的是二維模型,但材料參數和邊界條件是一維的。通過這個簡單的模型,驗證了自主開發的熱彈性有限元程序的正確性。月表附近的月巖跟月壤的熱彈性參數有很大的區別,因此,二維模型,甚至三維模型是必要的,在月球觀測和實驗的最新數據的基礎上,研究月壤與月巖在橫向上和縱向上的差異性、月表地形的非均勻起伏和太陽光照吸收熱量的時空變化不均勻性的影響,這也是未來進一步探究月表附近的熱月震的重要方向之一。

5 結論

本文利用熱彈性耦合有限元并行程序對月表月壤特征厚度差異造成的溫度、位移和熱應力時空演化進行建模計算,得出以下結論：

1)當月壤特征厚度為2 cm時,月表及50 cm以淺的月壤溫度、垂直位移和水平正應力隨時間成周期性變化,赤道地區溫度在100～385 K之間變化,垂直位移在-0.017～0.012 m之間升降,水平正應力在-1.731～2.066 MPa之間壓縮與拉張變化。這些周期性變化的外因在于月表太陽吸收熱量Qs和月球表面向外輻射的熱量QT相疊加的熱流邊界條件,同時月壤隨深度、溫度變化的熱力學參數的非線性是其內因。

2)下午到18：00日落時分,是月表溫度下降最快的時刻,同時也是月表水平正應力從壓應力轉為張應力變化速率最快增加的時刻。06：00日出時分,是月表溫度最低的時刻,同時也是月表水平拉應力最大的時刻。06：00日出后的時段,是月表溫度增加最迅速的時段,也是月表水平應力從最大拉張應力迅速轉變為壓應力的時段。日落和日出的相近時間,可能是月表附近熱事件或熱月震相對集中的時刻。更精確的模擬,還應考慮月表地形的起伏和月表附近介質的橫向不均勻性。

3)月壤溫度、垂直位移和水平正應力在月表以下較淺的深度變化較大。與零特征厚度模型相比,月壤特征厚度越大,無論是溫度、水平正應力還是垂直位移,差異也越大,并且不同時段的差異還有所不同,表明月壤特征厚度的影響需要進一步的定量研究。