標(biāo)準(zhǔn)線性固體模型的穩(wěn)定性條件數(shù)值解的計算方法及系統(tǒng)的制作方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及地震勘探基礎(chǔ)應(yīng)用技術(shù)領(lǐng)域,尤其涉及一種標(biāo)準(zhǔn)線性固體模型的有限 差分解的穩(wěn)定性條件數(shù)值解的計算方法及系統(tǒng)。
【背景技術(shù)】
[0002] 地震數(shù)值模擬是地震勘探和地震學(xué)的重要基礎(chǔ),同時也是了解復(fù)雜介質(zhì)中地震波 傳播規(guī)律的重要工具,其作用貫穿于整個地震采集、處理和解釋中。隨著地震勘探開發(fā)的深 入,常規(guī)的彈性介質(zhì)理論難W滿足實際介質(zhì)需求。地震波在實際地層中傳播時,能量和相位 都發(fā)生改變,直接影響地震資料的分辨率,實際地層表現(xiàn)出一定的黏彈特性,因此通過對黏 彈介質(zhì)中的地震波進行數(shù)值模擬,來研究和分析地震波傳播過程中的衰減特征,對實際地 震資料分辨率的提高非常有意義。但是確保黏彈介質(zhì)數(shù)值模擬穩(wěn)定(即時間步長的確定) 是一個急需解決的問題,送關(guān)系到模擬算法的成敗,也關(guān)系到模擬效率的高低。
[0003] 目前主要運用黏彈固體模型來表征實際介質(zhì)的粘滯特性,應(yīng)用最廣的黏彈固體模 型主要包括;Kelvin-Voigt固體模型、Maxwell固體模型和標(biāo)準(zhǔn)線性固體模型(也稱為標(biāo)準(zhǔn) 線性黏彈固體模型)。其中,Kelvin-Voigt固體模型不能考慮應(yīng)力作用下應(yīng)變的突然變化, 也不能表示應(yīng)力消失后的剩余應(yīng)變,Maxwell固體模型不具備蠕變特征,兩者都不足W描述 大多數(shù)黏彈介質(zhì)的特征,但是標(biāo)準(zhǔn)線性固體模型可W同時考慮具備應(yīng)變突然變化和剩余應(yīng) 變及蠕動特征,因此,標(biāo)準(zhǔn)線性固體模型能更加全面表征實際地層的黏彈特性,相比較而言 更符合實際情況,所W進行標(biāo)準(zhǔn)線性固體模型的黏彈介質(zhì)數(shù)值模擬是非常有必要的。
[0004] 現(xiàn)有技術(shù)中標(biāo)準(zhǔn)線性固體模型的黏彈介質(zhì)數(shù)值模擬主要是全Η維、全波場的模 擬,此種方法的計算量非常大,從而產(chǎn)生的費用(主要體現(xiàn)在電力、存儲設(shè)備、人力等資源 的消耗)是巨大的,數(shù)值模擬的效率低,因此有必要提供一種計算方法定量化地給出標(biāo)準(zhǔn) 線性固體模型數(shù)值模擬時間步長(即標(biāo)準(zhǔn)線性固體模型的有限差分解的穩(wěn)定性條件數(shù)值 解),從而更高效、更成功地指導(dǎo)數(shù)值模擬的進行。
【發(fā)明內(nèi)容】
[0005] 本發(fā)明所要解決的技術(shù)問題是現(xiàn)有技術(shù)中采用全Η維、全波場的模擬方法產(chǎn)生的 費用高、數(shù)值模擬效率低的缺陷。
[0006] 為了解決上述技術(shù)問題,本發(fā)明提供了一種標(biāo)準(zhǔn)線性固體模型的穩(wěn)定性條件數(shù)值 解的計算方法及系統(tǒng)。
[0007] 本發(fā)明的技術(shù)方案為:
[0008] -種標(biāo)準(zhǔn)線性固體模型的穩(wěn)定性條件數(shù)值解的計算方法,包括:
[0009] 構(gòu)建標(biāo)準(zhǔn)線性固體模型,并使所述標(biāo)準(zhǔn)線性固體模型包括彼此串聯(lián)的第一彈性體 和第二彈性體、W及與所述第一彈性體并聯(lián)的阻尼器;
[0010] 確定所述標(biāo)準(zhǔn)線性固體模型有限差分解的穩(wěn)定性條件的狀態(tài)傳遞矩陣;
[0011] 獲取多組參數(shù)集,并使每組參數(shù)集中包括所述第一彈性體的彈性系數(shù)、所述第二 彈性體的彈性系數(shù)、所述阻尼器的黏滯系數(shù)、頻率和介質(zhì)密度;
[0012] 在設(shè)定的空間差分精度和空間網(wǎng)格步長下,依次利用各組所述參數(shù)集并通過計算 機計算使得所述狀態(tài)傳遞矩陣的特征值的模小于1的時間步長;
[0013] 確定計算得出的時間步長為標(biāo)準(zhǔn)線性固體模型的穩(wěn)定性條件數(shù)值解。
[0014] 優(yōu)選的是,所述確定所述標(biāo)準(zhǔn)線性固體模型有限差分解的穩(wěn)定性條件的狀態(tài)傳遞 矩陣包括:
[0015] 確定所述標(biāo)準(zhǔn)線性固體模型的本構(gòu)方程
其中P為總應(yīng)力,ε為總應(yīng)變,彈性系數(shù)為Ml的 所述第一彈性體的應(yīng)變與黏滯系數(shù)為M2的所述阻尼器的應(yīng)變相等,Ms為所述第二彈性體的 彈性系數(shù);
[0016] 根據(jù)總應(yīng)變ε與質(zhì)點位移(u,v,w)間的關(guān)系方程
W及所述本 構(gòu)方程,得到第一方程:
[0017]
[0018] 對所述第一方程的左右兩邊分別對時間求二次偏導(dǎo)數(shù),得到第二方程:
[0019]
[0020] 利用所述第二方程和聲波的納維爾方程
得到第Η方程:
?
[0021] 其中Ρ為所述介質(zhì)密 度;
[0022] 將所述第S方程中的應(yīng)力取空間傅里葉變換,得到第四方程:
其中#為總應(yīng)力Ρ的空間傅里葉變換, k為波數(shù);
[0023] 對所述第四方程中的時間偏導(dǎo)數(shù)用差分近似,得到第五方程:
[0024]
其 中護聲"、;ri分別為第n-2, n-1,η, n+1時刻的夢值,并且所述波數(shù)k滿足:
在所述空間差分精度為2N的情況下,X,y,Z Η個方向上的空 間網(wǎng)格步長分別為Δχ,Ay, Δζ,曰1為對應(yīng)所述空間差分精度2Ν的空間差分系數(shù),At為 所述時間步長;
[0025] 根據(jù)所述第五方程,得到所述標(biāo)準(zhǔn)線性固體模型有限差分解的穩(wěn)定性條件的狀態(tài) 傳遞矩陣
其中:
[0029] 優(yōu)選的是,所述方法還包括;在計算得出的所述時間步長后,根據(jù)所述參數(shù)集計算 所述標(biāo)準(zhǔn)線性固體模型的品質(zhì)因子,并根據(jù)所述參數(shù)集和所述品質(zhì)因子計算所述標(biāo)準(zhǔn)線性 固體模型的介質(zhì)速度。
[0030] 優(yōu)選的是,在設(shè)定的空間差分精度和空間網(wǎng)格步長下,依次利用各組所述參數(shù)集, 運用安裝在所述計算機上的Matlab仿真軟件編程計算使得所述狀態(tài)傳遞矩陣的特征值的 模小于1的時間步長。
[0031] 一種標(biāo)準(zhǔn)線性固體模型的穩(wěn)定性條件數(shù)值解的計算系統(tǒng),包括:
[0032] 模型構(gòu)建單元,用于構(gòu)建標(biāo)準(zhǔn)線性固體模型,并使所述標(biāo)準(zhǔn)線性固體模型包括彼 此串聯(lián)的第一彈性體和第二彈性體、W及與所述第一彈性體并聯(lián)的阻尼器;
[0033] 狀態(tài)傳遞矩陣確定單元,用于確定所述標(biāo)準(zhǔn)線性固體模型有限差分解的穩(wěn)定性條 件的狀態(tài)傳遞矩陣;
[0034] 參數(shù)集獲取單元,用于獲取多組參數(shù)集,并使每組參數(shù)集中包括所述第一彈性體 的彈性系數(shù)、所述第二彈性體的彈性系數(shù)、所述阻尼器的黏滯系數(shù)、頻率和介質(zhì)密度;
[0035] 時間步長計算單元,用于在設(shè)定的空間差分精度和空間網(wǎng)格步長下,依次利用各 組所述參數(shù)集并通過計算機計算使得所述狀態(tài)傳遞矩陣的特征值的模小于1的時間步長;
[0036] 穩(wěn)定性條件數(shù)值解確定單元,用于確定計算得出的時間步長為標(biāo)準(zhǔn)線性固體模型 的穩(wěn)定性條件數(shù)值解。
[0037] 優(yōu)選的是,所述狀態(tài)傳遞矩陣確定單元包括:
[003引本構(gòu)方程確定單元,用于確定所述標(biāo)準(zhǔn)線性固體模型的本構(gòu)方程
其中P為總應(yīng)力,ε為總應(yīng)變,彈性系數(shù)為Ml的 所述第一彈性體的應(yīng)變與黏滯系數(shù)為M2的所述阻尼器的應(yīng)變相等,Ms為所述第二彈性體的 彈性系數(shù);
[0039] 第一方程確定單元,用于根據(jù)總應(yīng)變ε與質(zhì)點位移(u,v,w)間的關(guān)系方程
W及所述本構(gòu)方程確定單元確定的本構(gòu)方程,得到第一方程:
[0040]
[0041] 第二方程確定單元,用于對所述第一方程確定單元得到的第一方程的左右兩邊分 別對時間求二次偏導(dǎo)數(shù),得到第二方程:
[0042]
[0043] 第Η方程確定單元,用于利用所述第二方程確定單元得到的第二方程和聲波的納 維爾方程
得到第Η方程:
>
[0044] 其中Ρ為所述介質(zhì)密 度;
[0045] 第四方程確定單元,用于將所述第Η方程確定單元得到的第Η方程中的應(yīng)力取空 間傅里葉變換,得到第四方程:
其中# 為總應(yīng)力Ρ的空間傅里葉變換,k為波數(shù);
[0046] 第五方程確定單元,用于對所述第四方程確定單元得到的第四方程中的時間偏導(dǎo) 數(shù)用差分近似,得到第五方程:
[0047]
其中r-氣護1.護擴'' 分別為第η-化-lnn+1時刻的#值,并且所述波數(shù)k滿足:
在所述空間差分精度為2N的情況下,X,y,Z Η個方向上的空 間網(wǎng)格步長分別為Δχ,Ay, Δζ,曰1為對應(yīng)所述空間差分精度2Ν的空間差分系數(shù),At為 所述時間步長;
[0048] 狀態(tài)傳遞矩陣子確定單元,用于根據(jù)所述第五方程確定單元得到的第五方程,得 到所述標(biāo)準(zhǔn)線性固體模型有限差分解的穩(wěn)定性條件的狀態(tài)傳遞矩陣
癢中:
[0052] 優(yōu)選的是,所述系統(tǒng)還包括:
[0053] 品質(zhì)因子計算單元,用于在所述時