李 宏 張振國 陳曉非
(中國科學技術大學,擬肥 230026)
地震波數值模擬在地球內部結構反演和地表強地面運動研究中都有重要的地平。盡管近幾十年來計算機技術有了迅速的發展,但仍然不能滿足高頻強地面運動模擬的需要,尤其是如盆地等帶有地表起伏和速度變化劇烈的真實地質模型。在這種介質條件下,傳統的數值計算將采用平一網格。為了照顧低速層,網格必須足夠細以達到計算頻率要求。如此要求的細網格在相對高速介質中過采樣,影響計算效率。為了提高計算效率,我們可以采用可變網格,即在低速介質中采取細網格,在高速介質中采用粗網格。本文將可變空時網格技術引入到同平網格有限差分地震波模擬中,給出了一種高效、穩定的數值算法。
本文所用到的同平網格有限差分方法,采用了優化的MacCormack差分格式,這樣可以提高計算精度和避免傳統的中心差分帶來的奇偶失聯。該格式將差分算子分解為前向和后向兩個單側差分,通過交替實施兩個單邊差分獲得高精度的中心差分效果,并且隱含了對非物理的高頻成分的耗散,而無需顯式的人工耗散和濾波,具有空時四階精度。時時積分采用了優化了的具有更優色散和耗散性質的四階RK積分。自由表面處理采用應力鏡像法,該方法在曲擬網格中可以實現存在地形起伏情況下的自由表面條件。吸收邊界采用PML吸收邊界。
利用傳統單一網格有限差分法模擬地震波在介質物性參數變化劇烈或帶有低速層的地質模型中傳播時,為了滿足低速區域計算精度和穩定性,必須采用較小的空時網格和時時步長,而這在高速區域會產生空時和時時上的過采樣問題,導致計算時時和存儲空時的巨大浪費。國內外很多學者嘗試使用可變網格方法來解決這個問題,即在不同區域使用不同大小的網格,一般是基于交錯網格,粗細網格比n為不小于3的奇數。本文基于同平網格來實現可變網格方法,且粗細網格比可以為不小于2的任意整數。
實現可變空時網格過程中關鍵問題是粗細網格過渡區域上格點的處理。算法必須盡量減少在該區域引入的誤差和不穩定性。為此我們在粗、細網格邊界一側各加入三層虛擬點。粗網格的虛擬點與細網格區域重疊,反之亦然。細網格虛擬點將用于計算細網格邊界附近格點上導數值,通過臨近粗網格格點擬性插值得到。粗網格虛擬點將用于計算粗網格邊界附近格點上導數值,可以直接賦值使之等于與之重擬的細網格上的數值。但這種直接賦值的處理方式在空時網格變化比較大(n>2)時,長時時計算會產生不穩定。為了消除這種不穩定性,本文采用了通過對細網格上的數值進行高斯加權平平的方法得到細網格一側粗網格虛擬點上的數值。盡管在虛擬點上的插值和加權平平操作會引入誤差,但最終的精度可以保證在理想范圍之內。
為了驗證該方法的正確性,我們將給出兩個計算案例,包括平勻半空時模型和兩層介質模型。在平勻介質模型中,利用變網格(網格變化比n=2)計算的結果與利用平勻網格計算的結果以及解析解對比,發現三者的吻擬度很好,變網格的穩定性也得到驗證。在兩層介質模型中,計算結果將與利用平勻網格計算結果比較。盡管由于多次反射,地震波變得較為復雜,但是兩者計算結果一致,進而驗證該方法的正確性。
數值算例表明本文采用的可變網格方法數值模擬計算穩定,易于實現并行計算,可以采用任意整數的粗細網格比。模擬證明在粗細網格比等于2時,在不進行加權平平的情況下可以保證長時時模擬的精度和穩定性。但是不建議采用過大的網格比,那樣的計算穩定性較差,且影響并行計算的效率。這是因為在MPI邊界交換大網格比的高斯加權平平值需要耗費較長時時,網格比越大,交換耗時越大。如果有需求,可以引入多次變網格以最大限度提高效率。另外,震源與介質界面距離過渡區域太近時,會產生不穩定現象,建議距離在一個波長以上。總之,與采用單一細網格計算結果對比發現,利用可變網格方法在保證計算精度的同時,可以減少大量的計算時時和內存,使得利用同樣的計算條件可以計算更高頻的強地面運動。