[發明專利]復雜介質彈性波傳播模擬的多網格切比雪夫并行譜元法

有效

• 申請號: 201310452393.5
• 申請日: 2013-09-27
• 公開（公告）號: CN103530451B
• 公開（公告）日: 2016-09-28
• 發明（設計）人: 蘇暢;G·塞里安尼;林偉軍
• 申請（專利權）人: 中國科學院聲學研究所
• 主分類號: G06F17/50
• 分類號: G06F17/50
• 代理公司: 北京億騰知識產權代理事務所 11309
• 代理人: 陳霽

技術領域

本發明涉及地球物理，工程地震學，計算聲學等領域的高性能計算。

背景技術

復雜介質中彈性波的傳播模擬在地球物理，工程地震學和計算聲學等領域具有重要地位。如何找到一種更精確、更快速、更適宜并行化的數值模擬方法也一直是國、內外相關領域研究者的工作目標和重點。

當前在彈性波方程模擬中普遍使用的數值方法，主要有有限差分法、有限元法、偽譜法和譜元法等。它們有各自的優缺點并在不同的領域得到應用。

有限差分法將方程中的導數用差分近似，從而將偏微分方程轉換為代數方程來求解，其誤差由離散化配置點數目和Taylor級數的截斷誤差決定。有限差分法由于編程簡捷和計算速度較快被廣泛應用于計算地球物理和工程地震中。但低階有限差分法需要大量的網格點來保證精度，而高階有限差分法又很難高效地處理自由界面和復雜結構問題。

有限元方法基于波動方程的弱形式，將計算區域劃分為有限個互不重疊的單元，在每個單元內，波場由低階多項式(如分段線性函數)來近似。有限元法適合處理邊界和內部復雜結構，在結構分析和瞬態模擬方面獲得了廣泛的應用。但有限元法在大型的彈性波方程模擬領域應用相對較少，其原因一是有限元法需要耗用大量計算資源，二是低階有限元存在頻散現象，而通常的高階有限元存在偽波。

偽譜法利用無限可微的全局函數（如富里葉Fourier級數），將問題從物理空間轉化為波數域進行處理，具有良好的處斂性。偽譜法只需要很少的空間網格點數即可達到很高的精度。為了處理邊界條件，在空間展開時經常使用切比雪夫Chebyshev或拉格朗日Legendre正交多項式來取代Fourier級數。但由于正交多項式的配置點分布很不均勻，給時間步長的選擇帶來很大限制。同時，如同有限差分法一樣，偽譜法不能自如地處理彎曲的自由界面和復雜結構的問題。為解決這類問題有人提出了彎曲坐標和區域分解等方法，但代價是計算量的增加。

譜元法由Patera在1984年提出，早期主要應用于流體力學。Seriani等在1991年首次將Chebychev譜元法引入到聲波方程的數值模擬中，各國研究者在其后做了大量的相關研究。它把有限元法和偽譜法相結合，兼具了有限元處理邊界和結構的靈活性和偽譜法的快速收斂特性。在達到相同的精度前提下，可以采用比傳統有限元更稀疏的單元劃分，減少了計算時間和內存需求。

譜元法的基本方法是在每一個單元上用高階譜展開。選取以截斷的正交多項式表示的基函數，在各個單元上利用配置點插值，以提高解的收斂速度。其主要步驟是：(1)首先把計算的區域分成許多子域（單元），每個子域由若干節點（配置點）組成;(2)在每個子域中把近似解表示成截斷的正交多項式展開;(3)用Galerkin方法求解正交問題的變分格式，得到全局的近似解。

譜元法可以用較稀疏的網格和單元獲得較高精度，但傳統的譜元法中，每個單元內只能有單一均勻介質，某些情況下嚴重降低計算效率。比如當介質結構復雜，變化尺度很小甚至小于波長時，必須按介質變化的小尺度劃分單元和求解，這樣將會造成極大的計算規模。

此外，傳統的譜元法求解波動方程時，需要形成全局剛度矩陣和質量矩陣，考慮衰減時還需引進阻尼矩陣，需要耗費大量的存貯空間，在算法上限制了其并行化的效率。注意到在波動方程的時間迭代過程中，需要做多次矩陣和向量乘積，如Sp_k，這里S和p_k分別為全局剛度（或質量）矩陣和全局向量。考慮到全局矩陣的稀疏特性，全局的矩陣和相量相乘耗費了大量的計算時間。

以上兩點都極大地阻礙了譜元法在地球物理，工程地震學，計算聲學等領域的大規模實際應用。

發明內容

本發明的目的在于提供一種能夠克服上述缺點的算法。

為此，本發明提供一種多網格Chebyshev并行譜元方法。該方法包括：

a.將計算區域劃分為大而規則的單元；

b.在單元內定義主網格和輔助網格，主網格為Chebyshev配置點；

c.利用主網格將單元內波場作Chebyshev截斷展開，利用輔助網格將單元內介質參數和外力作截斷展開；

d.將波場、介質參數和外力的截斷展開式代入波動方程，得到單元質量矩陣、單元剛度矩陣和單元外力向量；

e.基于各單元的單元質量矩陣、單元剛度矩陣和單元外力向量，在主網格上求解彈性波方程。

優選地，利用預條件共軛梯度法求解彈性波方程。