飽和兩相介質(zhì)近場波動問題的一種時域全顯式數(shù)值計算方法

李 亮，李 果，杜修力，宋 佳

(1. 北京工業(yè)大學(xué)城市與工程安全減災(zāi)教育部重點實驗室，北京 100124；2. 北方工業(yè)大學(xué)土木工程學(xué)院，北京 100144)

在沿海、沿江地區(qū)的各類工程場地中廣泛存在著飽和土地基，許多工程結(jié)構(gòu)，如跨海、越江隧道以及地鐵車站等地下結(jié)構(gòu)也修建于飽和土地基之中。飽和土可抽象為飽和兩相介質(zhì)進(jìn)行相關(guān)力學(xué)問題的研究，這是一種由固相(土顆粒構(gòu)成的骨架)與液相(土骨架孔隙中的流體)共同構(gòu)成的二元混合介質(zhì)。土動力學(xué)與巖土地震工程中的一些典型問題，諸如地震波在飽和土場地中的傳播，結(jié)構(gòu)與飽和土的動力相互作用等問題均可歸結(jié)為飽和兩相介質(zhì)的近場波動問題進(jìn)行研究。進(jìn)行近場波動問題的研究時，需將求解區(qū)域劃分為兩部分，結(jié)構(gòu)及其鄰近的場地為內(nèi)域，內(nèi)域以外的場地為外域。飽和兩相介質(zhì)近場波動問題的時域求解需解決兩方面問題[1]：一是內(nèi)域數(shù)值算法的建立，即建立內(nèi)域中飽和兩相介質(zhì)動力響應(yīng)的時域計算方法；二是外域?qū)?nèi)域的動力作用的描述方法，即所謂的人工邊界條件建立。

飽和兩相介質(zhì)的近場波動問題是由描述波動規(guī)律的波動方程與相應(yīng)的定解條件(邊界條件和初始條件)共同構(gòu)成的定解問題。建立該問題的內(nèi)域數(shù)值計算方法需針對特定形式的波動方程進(jìn)行。根據(jù)方程中包含的基本未知量的不同，飽和兩相介質(zhì)的波動方程包括：1)以固相位移和液相位移為基本未知量的u-U 形式方程[2]；2)以固相位移和固相與液相的相對位移為基本未知量的u-w 形式方程[3]；3)以固相位移、固相與液相的相對位移以及液相(孔隙流體)壓力為基本未知量的u-w-p 形式方程[4]；4)以固相位移和孔隙流體壓力為基本未知量的u-p 形式方程[5]。該方程是在u-U 形式方程的基礎(chǔ)上略去了固相與液相的相對加速度項，即忽略兩相之間的慣性的差異而得到的簡化形式的波動方程。該方程適用于描述包含地震響應(yīng)在內(nèi)的低頻動力問題。作為一種由矢量未知量和標(biāo)量未知量共同構(gòu)成的混合形式的方程，u-p 形式波動方程相對于其他全矢量形式的方程具有相對較少的未知量數(shù)目，可以在一定程度上減少計算求解的工作量；基于該方程可以直接求解得到孔隙流體壓力p，其是很多土動力學(xué)與巖土地震工程問題(如液化判別問題)中需要求解的關(guān)鍵變量之一?；谏鲜鰞煞矫娴脑颍瑄-p 形式的波動方程在飽和兩相介質(zhì)的近場波動問題中得到了廣泛的應(yīng)用。

由于方程形式的特殊性，針對u-p 形式波動方程的時域數(shù)值計算方法的基本模式為交替求解，即在一個時間分步上采用不同的計算模式分別計算固相位移u 和孔隙流體壓力p，且計算過程須按照一定的順序進(jìn)行。根據(jù)求解模式的特點，交替解法可以分為隱式方法、顯式方法和二者相結(jié)合的顯-隱式方法。隱式方法的特點是采用求解耦聯(lián)方程組的模式計算固相位移u 和孔隙流體壓力p，獲得飽和兩相介質(zhì)動力響應(yīng)的時域數(shù)值結(jié)果；而顯式方法的特點是無須求解耦聯(lián)方程組，采用迭代或遞推的計算模式完成動力響應(yīng)時程的計算，獲得飽和兩相介質(zhì)動力響應(yīng)的時域數(shù)值結(jié)果。Park[6]、Zienkiewicz 等[7]、黃茂松等[8 ? 9]和李曉欣等[10]分別建立了針對u-p 形式波動方程的隱-隱式交替的時域逐步積分求解算法，這些算法對固相位移u 和孔隙流體壓力p 均采用隱式的求解模式進(jìn)行時域的計算。隱式方法具有良好的計算穩(wěn)定性，但是由于在計算的每一個時間分步上都需要求解耦聯(lián)方程組，對于自由度數(shù)目巨大的動力問題其計算效率將存在一定的問題。顯-隱式交替解法也可稱為半顯式解法，是對隱式方法的改進(jìn)與發(fā)展，其特點是顯式與隱式兩種時域計算模式的結(jié)合，即針對需要求解的兩個動力響應(yīng)中的某一個，如固相位移u 采用顯式計算模式，而另一個動力響應(yīng)的求解則采用隱式計算模式。由于顯式計算模式的部分引入，可以使計算效率得到一定程度的改善與提高。Zienkiewicz 等[11 ? 12]、Pastor 等[13]、Li 等[14]和杜修力等[15]分別建立了基于u-p 形式波動方程的飽和兩相介質(zhì)近場波動問題的顯-隱式時域交替解法。

對于飽和兩相介質(zhì)近場波動問題的時域求解，尤其是動力自由度數(shù)目巨大的復(fù)雜問題，全顯式的時域計算方法是較好的解決方案，即針對固相位移和孔隙流體壓力均采用迭代與遞推的求解模式。由于完全避免了求解耦聯(lián)方程組的過程，因此全顯式時域計算方法的求解效率相對于隱式方法和顯-隱式方法有比較顯著的改進(jìn)與提升。國內(nèi)學(xué)者針對飽和兩相介質(zhì)近場波動問題的時域全顯式計算方法開展了一定的研究工作。段雪銘等[16]將精細(xì)時程積分技術(shù)引入飽和兩相介質(zhì)動力問題的求解，建立了針對該問題的一種新的時域計算方法。由于引入了精細(xì)時程積分技術(shù)，該方法具有較好的計算精度；同時，該方法的求解過程以迭代的模式完成，避免了矩陣求逆，具有時域顯式計算方法的基本特點。但是該方法的計算過程中需進(jìn)行矩陣的指數(shù)運算，使其計算效率受到影響。在近場波動問題時域顯式計算方法的建立過程中，為得到解耦的時域顯式迭代計算格式，須對波動方程的系數(shù)矩陣進(jìn)行對角化處理以消除動力耦聯(lián)。宋佳等[17 ? 18]在波動方程系數(shù)矩陣對角化的基礎(chǔ)上，建立了兩種飽和兩相介質(zhì)近場波動問題的時域全顯式求解方法。第一種方法中分別應(yīng)用中心差分法和精細(xì)時程積分技術(shù)計算固相和液相的時域動力響應(yīng)，而第二種方法中分別基于杜修力和王進(jìn)廷[19]提出的時域顯式算法和單邊差分法計算固相位移和孔隙流體壓力。

本文將進(jìn)一步開展針對飽和兩相介質(zhì)近場波動問題的全顯式時域求解方法的研究工作，綜合應(yīng)用中心差分法、向后差分法以及Newmark 常平均加速度法，建立飽和兩相介質(zhì)近場波動問題的一種新的時域全顯式數(shù)值計算方法。將該方法的數(shù)值結(jié)果與相應(yīng)的解析解結(jié)果進(jìn)行對比，完成方法的驗證。并將該方法與透射人工邊界[20]相結(jié)合，應(yīng)用于飽和兩相介質(zhì)近場波動問題的研究之中。

1 u-p 形式的飽和兩相介質(zhì)波動方程

假定飽和兩相介質(zhì)中的固相為各向同性的線彈性介質(zhì)，其應(yīng)力-應(yīng)變關(guān)系服從胡克定律；液相(孔隙流體)為可壓縮的理想流體，并且在忽略體積力的情況下，以固相位移u 和孔隙流體壓力p 作為基本未知量的u-p 形式的飽和兩相介質(zhì)彈性波動方程可以表達(dá)為如下的形式：

式中， λ和G 分別為固相的Lame 常數(shù)和剪切模量。

梯度算子向量 ?和微分算子矩陣 L的表達(dá)式分別為：

向量 m的表達(dá)式為：

已有的研究工作表明[21]，式(2)中與固相加速度相關(guān)的第三項對計算結(jié)果的影響很小，故在后續(xù)的研究中略去該項。

應(yīng)用伽遼金法[22]對波動方程式(1)和式(2)進(jìn)行有限元空間離散，得到其伽遼金弱式如下：

式中： M 為質(zhì)量矩陣； K 為剛度矩陣；Q 為固相與孔隙流體動力響應(yīng)的耦合矩陣； J為孔隙流體的滲流矩陣， S 為孔隙流體的壓縮矩陣； fu和 fp分別為作用在固相與孔隙流體上的邊界荷載向量。式(7)和式(8)中各個矩陣和向量的表達(dá)式為：

式中： Nu和 Np分別為針對固相位移u 和孔隙流體壓力p 的插值形函數(shù)； fue和 fpe分別為作用在邊界單元上對應(yīng)于固相和孔隙流體的荷載向量。

2 方法的建立

2.1 波動方程系數(shù)矩陣對角化

經(jīng)有限元空間離散后得到的飽和兩相介質(zhì)波動方程式(7)和式(8)為具有動力耦聯(lián)的耦合方程組，其各個系數(shù)矩陣描述有限元計算模型中的各節(jié)點之間的動力相互作用，均為非對角矩陣。欲構(gòu)建針對該方程組的時域顯式求解方法，建立動力響應(yīng)的時域顯式迭代計算公式，須對該方程組進(jìn)行解耦，解耦的方法是對方程組的某些系數(shù)矩陣進(jìn)行對角化處理，得到對角形式的系數(shù)矩陣，從而可以消除不同節(jié)點動力響應(yīng)之間的耦聯(lián)。對角化后的系數(shù)矩陣的求逆運算只須針對主對角元素進(jìn)行，避免了求解過程中的矩陣整體求逆，而這正是影響計算效率的主要因素之一。

在本文方法中，需對質(zhì)量矩陣 M和孔隙流體壓縮矩陣 S進(jìn)行對角化處理，得到集中質(zhì)量矩陣和對角形式的孔隙流體壓縮矩陣。基于對有限元計算系統(tǒng)動力特性的認(rèn)識引入相關(guān)的假定：1)由于有限元計算系統(tǒng)中各節(jié)點之間的慣性耦合相對于節(jié)點慣性本身為高階小量，因此可以忽略各節(jié)點之間的慣性耦合，認(rèn)為各節(jié)點的慣性相互獨立，將單元質(zhì)量集中到各個節(jié)點上；2)忽略有限元計算系統(tǒng)中各節(jié)點的流體壓力時間變化率之間的耦合，認(rèn)為各節(jié)點的流體壓力時間變化率相互獨立，與流體壓力時間變化率相關(guān)的矩陣 S的非對角元素為零?；谏鲜黾俣ǎo出矩陣對角化的方法如下[22]：

對角化后的單元質(zhì)量矩陣的主對角元素為各節(jié)點上的集中質(zhì)量，各節(jié)點上的集中質(zhì)量之和等于單元的總質(zhì)量；非對角元素為零，即形成集中質(zhì)量矩陣。對角化后的單元質(zhì)量矩陣的表達(dá)式為：

應(yīng)用相同的方法可完成流體壓縮矩陣 S 的對角化，對角化后的單元流體壓縮矩陣的表達(dá)式為：

式(17)和式(19)中： Ve為單元的體積； Id為單位向量；d 為單元在幾何空間的維數(shù)。

2.2 動力響應(yīng)的時域顯式逐步計算列式

基于系數(shù)矩陣對角化后的飽和兩相介質(zhì)波動方程組，推導(dǎo)建立飽和兩相介質(zhì)動力響應(yīng)的時域顯式逐步積分計算的計算列式。將動力荷載(如地震荷載)的作用時間t 按照固定的步長 ?t進(jìn)行離散，得到一系列離散的時間點，第k個時間點對應(yīng)的時刻為tk=k?t, (k=1,2,3···)，k 也是迭代計算的時間步數(shù)。本文方法中動力響應(yīng)時域計算的基本模式為由第k步(tk時刻)的動力響應(yīng)遞推得到第k+1步(tk+1時刻)的動力響應(yīng)。

迭代計算的時間步數(shù)為k時(對應(yīng)的計算時刻為tk)的波動方程組可以表示為如下的形式：

式中， Ml和 Sl分別為對角化后的整體質(zhì)量矩陣和流體壓縮矩陣。

基于中心差分法，第k 步的固相加速度的計算式為：

將式(22)代入式(20)，得到第k+1步的固相位移的表達(dá)式如下：

假定在時間步長? t的范圍內(nèi)，流體壓力 p為 線性變化，其對時間的變化率(導(dǎo)數(shù))可以用有限差分的形式近似表示。采用向后差分法，即：

將式(24)代入式(21)，得到第k+1步的孔隙流體壓力的計算式如下：

迭代計算的時間步數(shù)為k+1時(對應(yīng)的計算時刻為tk+1)的式(20)可以寫成如下形式：

將式(20)與式(26)進(jìn)行疊加，可得：

Newmark 常平均加速度法采用如下假定：

將式(28)代入式(27)，得到第k+1步的固相速度的表達(dá)式如下：

式(23)、式(25)和式(29)共同組成了本文算法中進(jìn)行動力響應(yīng)計算求解的時域顯式逐步積分列式。應(yīng)用該組計算列式，可以按照計算時間步數(shù)依次進(jìn)行迭代求解，完成飽和兩相介質(zhì)動力響應(yīng)時程的計算。由于質(zhì)量矩陣 Ml和流體壓縮矩陣Sl均為對角矩陣，上述計算列式中的矩陣求逆只須除以主對角線上的元素，無須進(jìn)行矩陣的整體求逆。本文算法的求解過程通過迭代遞推的方式完成，具有時域顯式數(shù)值計算方法的基本特征。

在迭代計算的每個時間分步上，動力響應(yīng)的求解須按照如下步驟進(jìn)行：1)分別應(yīng)用式(23)和式(25)由前一個時間分步的動力響應(yīng)計算當(dāng)前分步的固相位移uk+1和孔隙流體壓力pk+1；2)應(yīng)用式(29)由前一時間分步的動力響應(yīng)以及步驟1)得到的uk+1和 pk+1計算當(dāng)前分步的固相速度u˙k+1。按照計算時間步數(shù)對上述求解步驟進(jìn)行循環(huán)，即可完成飽和兩相介質(zhì)動力響應(yīng)時程的計算?；诒疚乃惴ǖ那蠼饽Ｊ娇芍?，飽和兩相介質(zhì)波動問題相關(guān)的全部動力響應(yīng)，即固相位移 u 和固相速度 u˙，以及孔隙流體壓力 p都是采用顯式的計算模式，即迭代和遞推的模式進(jìn)行時域求解，因此本文建立的算法是一種全顯式的時域數(shù)值求解方法。

3 方法的驗證

3.1 矩陣對角化合理性驗證

如前文所述，本文方法的建立過程中，為了得到動力響應(yīng)的時域顯式迭代計算公式，將質(zhì)量矩陣和流體壓縮矩陣進(jìn)行了對角化處理，將原有的非對角形式的一致質(zhì)量矩陣和流體壓縮矩陣轉(zhuǎn)化為對角形式的集中質(zhì)量矩陣和流體壓縮矩陣。本小節(jié)針對具體的算例進(jìn)行矩陣對角化合理性的驗證。算例的計算模型為飽和土場地，場地的長度和深度均取30 m。采用四節(jié)點矩形單元對計算區(qū)域進(jìn)行有限元空間離散，單元尺寸為3 m×3 m。計算模型的上表面設(shè)置為自由排水邊界；左、右兩側(cè)設(shè)置為不透水邊界，并約束其水平位移；模型底部邊界設(shè)置為不透水的固定邊界。計算模型如圖1 所示。在計算模型的上表面施加均布的正弦荷載q=1+0.5sin(2t) kPa，持續(xù)時長為30 s，荷載時程如圖2 所示。

圖1 飽和土場地動力響應(yīng)計算模型Fig.1 Calculation model for dynamic response of saturated soil field

圖2 正弦荷載時程Fig.2 Time history of applied sinusoidal load

將坐標(biāo)原點設(shè)置在計算模型的左下角，選取模型中的三個節(jié)點A(3, 3)、B(3, 15)、C(3, 27)計算其動力響應(yīng)。飽和土場地的材料參數(shù)列于表1中，計算的時間步長取0.0001 s。

表1 飽和土場地材料參數(shù)Table1 Material parameters of saturated soil field

文獻(xiàn)[23]中提出了針對u-p 形式波動方程的飽和兩相介質(zhì)動力問題的顯-隱式計算方法，文中算例表明該方法具有良好的計算精度。將本文建立的基于矩陣對角化的時域顯式算法的計算結(jié)果與文獻(xiàn)[23]中方法的計算結(jié)果進(jìn)行對比，以驗證矩陣對角化的合理性和可行性。應(yīng)用本文算法得到的三個計算節(jié)點的固相位移時程與文獻(xiàn)[23]方法的相應(yīng)計算結(jié)果的對比如圖3 所示，三個計算節(jié)點的孔隙流體壓力時程的對比如圖4～圖6 所示。

圖3 計算節(jié)點固相位移時程對比Fig.3 Comparison of solid-phase displacement time history of calculated nodes

圖4 節(jié)點A 孔隙流體壓力時程對比Fig.4 Comparison of pore-pressure time historyof node A

圖5 節(jié)點B 孔隙流體壓力時程對比Fig.5 Comparison of pore-pressure time historyof node B

圖6 節(jié)點C 孔隙流體壓力時程對比Fig.6 Comparison of pore-pressure time historyof node C

由圖3～圖6 可知，本文基于矩陣對角化的時域顯式算法的計算結(jié)果與文獻(xiàn)[23]中方法的計算結(jié)果符合良好，這說明矩陣對角化在時域顯式數(shù)值計算方法的建立中具有合理性和可行性，基于矩陣對角化能夠建立針對飽和兩相介質(zhì)近場波動問題的正確的時域數(shù)值計算方法。

3.2 方法的解析解驗證

在飽和兩相介質(zhì)近場波動問題數(shù)值計算方法的研究中，對建立的數(shù)值算法進(jìn)行驗證的通行做法是解析解驗證，即針對特定的算例將算法的數(shù)值結(jié)果與相應(yīng)的解析解結(jié)果進(jìn)行對比。本小節(jié)進(jìn)行本文建立的時域顯式數(shù)值計算方法的解析解驗證。

算例為一維飽和土柱在表面均布荷載作用下的動力響應(yīng)。土柱的寬度為3 m；高度為200 m。上表面設(shè)置為自由排水邊界，左、右兩側(cè)設(shè)置為不透水邊界，并約束其水平位移；底面設(shè)置為不透水的固定邊界。邊界條件的表達(dá)式如下：

在土柱頂面施加均布的動力荷載。本算例的計算模型如圖7 所示。作用的動力荷載形式為階躍荷載和正弦荷載，荷載時程分別如圖8 和圖9所示。飽和土柱的材料參數(shù)見表2，表2 中參數(shù)均為無量綱化的參數(shù)。

圖7 飽和土柱動力響應(yīng)計算模型Fig.7 Calculation model for dynamic response of saturated soil column

圖8 正弦荷載時程Fig.8 Time history of applied sinusoidal load

圖9 階躍荷載時程Fig.9 Time history of applied step load

表2 飽和土柱材料參數(shù)Table2 Material parameters of saturated soil column

圖10 正弦荷載作用下飽和土柱上表面固相位移時程對比Fig.10 Comparison of solid-phase displacement time history of top surface of saturated soil column under sinusoidal loading

圖11 階躍荷載作用下飽和土柱上表面固相位移時程對比Fig.11 Comparison of solid-phase displacement time history of top surface of saturated soil column under step loading

由圖10 和圖11 可知，應(yīng)用本文算法得到的飽和土柱動力響應(yīng)的數(shù)值解結(jié)果與相應(yīng)的解析解結(jié)果符合良好，表明本文算法能夠?qū)︼柡蛢上嘟橘|(zhì)的動力響應(yīng)進(jìn)行較為準(zhǔn)確的計算，驗證了本文算法的正確性。

4 方法的應(yīng)用

本節(jié)將本文建立的時域顯式數(shù)值計算方法應(yīng)用于飽和兩相介質(zhì)的近場波動問題之中，進(jìn)行飽和土場地地震響應(yīng)的計算研究。本項研究的意義表現(xiàn)在兩個方面：第一，場地的地震響應(yīng)特性是土動力學(xué)和巖土地震工程領(lǐng)域中的重要研究課題之一[25]；第二，場地地震響應(yīng)的計算是地震動輸入條件下場地土體與結(jié)構(gòu)動力相互作用問題研究的基礎(chǔ)和重要步驟。

建立如圖12 所示的飽和土場地地震響應(yīng)計算模型。場地的高度和寬度均為20 m，采用四節(jié)點矩形單元進(jìn)行計算區(qū)域的有限元離散，單元尺寸為2 m×2 m，計算模型共包含100 個單元，121 個節(jié)點。計算模型的上表面設(shè)置為自由排水邊界，左、右兩側(cè)邊界和底部邊界采用透射人工邊界方法[20]進(jìn)行處理。透射人工邊界屬于局部人工邊界條件中的位移型人工邊界，其基本原理是通過對外行波傳播規(guī)律的直接模擬建立人工邊界上位移的時空外推計算公式。具體的實現(xiàn)方法為將邊界節(jié)點當(dāng)前時刻的動力響應(yīng)表達(dá)為該節(jié)點前幾個時刻的動力響應(yīng)以及其鄰近的內(nèi)節(jié)點前幾個時刻動力響應(yīng)的線性組合。基于透射人工邊界方法的基本原理，得到邊界節(jié)點1 在 p時刻的固相位移的表達(dá)式為：

圖12 飽和土場地地震響應(yīng)計算模型Fig.12 Calculation model for seismic response of saturated soil field

采用1995 年日本阪神地震中的Kobe 地震記錄作為輸入地震動，該地震記錄的持時為20 s，其位移時程和速度時程分別如圖13 和圖14 所示。令地震動分別以壓縮波(P 波)和剪切波(S 波)的形式從計算模型的底部邊界垂直輸入。取計算模型上表面的中間節(jié)點A 和模型深度中點位置的中間節(jié)點B 為計算節(jié)點。計算的時間步長為0.0001 s，計算輸入的飽和土場地的材料參數(shù)列于表3。

圖13 Kobe 地震記錄位移時程 Fig.13 Displacement time history of Kobe earthquake record

圖14 Kobe 地震記錄速度時程Fig.14 Velocity time history of Kobe earthquake record

表3 飽和土場地材料參數(shù)Table3 Material parameters of saturated soil field

當(dāng)?shù)卣鸩ㄒ訮 波形式入射時，應(yīng)用本文方法得到的飽和土場地表面(A 點)的固相位移時程和固相速度時程(豎直方向)的計算結(jié)果分別如圖15和圖16 所示，場地中部(B 點)的孔隙流體壓力時程如圖17 所示；當(dāng)?shù)卣鸩ㄒ許 波形式入射時，飽和土場地表面(A 點)的固相位移時程和固相速度時程(水平方向)的計算結(jié)果分別如圖18 和圖19所示，場地中部(B 點)的孔隙流體壓力時程如圖20 所示。對比圖15、圖16 和圖18、圖19 所示的飽和土場地表面的動力響應(yīng)時程與圖13 和圖14所示的入射地震動的時程可知，二者的波形變化趨勢保持一致。同時，動力響應(yīng)的幅值相對于相應(yīng)的入射地震動的幅值有近似為2 倍的放大關(guān)系。由圖17 和圖20 可知，飽和土場地中部的孔隙流體壓力時程與入射地震動時程的波形變化趨勢也基本相符。以上計算結(jié)果符合彈性波動理論的規(guī)律[1]，表明本文建立的時域顯式數(shù)值計算方法可以應(yīng)用于飽和土場地地震響應(yīng)問題的計算研究，得到正確的計算結(jié)果。上述問題屬于典型的飽和兩相介質(zhì)近場波動問題，因此，本文方法對飽和兩相介質(zhì)近場波動問題具有較好的適用性。

圖15 P 波入射場地表面固相位移時程Fig.15 Solid-phase displacement time history of field surface for P-wave incidence

圖16 P 波入射場地表面固相速度時程Fig.16 Solid-phase velocity time history of field surface for P-wave incidence

圖17 P 波入射場地中部孔隙流體壓力時程Fig.17 Pore-pressure time history of central part of field for P-wave incidence

圖18 S 波入射場地表面固相位移時程Fig.18 Solid-phase displacement time history of field surface for S-wave incidence

圖19 S 波入射場地表面固相速度時程Fig.19 Solid-phase velocity time history of field surface for S-wave incidence

圖20 S 波入射場地中部孔隙流體壓力時程Fig.20 Pore-pressure time history of central part of field for S-wave incidence

5 方法穩(wěn)定性研究

時域顯式數(shù)值計算方法具有條件穩(wěn)定性。對其進(jìn)行穩(wěn)定性特性研究的基本方法為對該方法的時域迭代計算格式相應(yīng)的傳遞矩陣進(jìn)行研究，由傳遞矩陣的譜半徑不大于1 的限定條件獲得該方法的穩(wěn)定性條件。

本文建立的時域全顯式數(shù)值計算方法的傳遞矩陣為：

記ρ(T)為傳遞矩陣T 的譜半徑，ρ(T)=max(λi)，λi為傳遞矩陣的特征值。當(dāng)ρ(T)≤1，且傳遞矩陣T 的特征方程的重根滿足|λi|≤1時，則傳遞矩陣T 是穩(wěn) 定的，即相應(yīng)的數(shù)值積分算法是穩(wěn)定的。為簡化分析，將飽和兩相介質(zhì)系統(tǒng)的物理力學(xué)參數(shù)取定值，則其波動方程中的系數(shù)矩陣均為常量矩陣，此時傳遞矩陣T 的系數(shù)行列式為：

令

由行列式|λE?T|=0得到傳遞矩陣T 的特征多項式如下：

式中：

式(37)可以改寫為如下形式：

式中：

若式(38)存在負(fù)實部解，須滿足Hurwitz 判據(jù)[26]。根據(jù)該判據(jù)可知，系統(tǒng)穩(wěn)定的充要條件為：當(dāng)多項式的所有系數(shù)均為正數(shù)時，多項式所有系數(shù)的主行列式及順序主子式全部為正。由上述判據(jù)可得到本文建立的時域全顯式數(shù)值算法的穩(wěn)定性條件為：

以時間步長表示的穩(wěn)定性條件為：

6 結(jié)論

本文基于以固相位移u 和孔隙流體壓力p 作為基本未知量的u-p 形式的飽和兩相介質(zhì)彈性波動方程，開展了飽和兩相介質(zhì)近場波動問題時域顯式數(shù)值計算方法的研究。通過對波動方程中的質(zhì)量矩陣和孔隙流體壓縮矩陣進(jìn)行對角化處理，消除了方程中的動力耦聯(lián)，實現(xiàn)了波動方程的解耦?；诮怦畹挠邢拊▌臃匠探M，分別應(yīng)用中心差分法和Newmark 常平均加速度法求解固相位移和速度，基于向后差分法求解孔隙流體壓力，推導(dǎo)得到了飽和兩相介質(zhì)動力響應(yīng)的時域顯式逐步積分計算的計算列式，建立了飽和兩相介質(zhì)近場波動問題的一種新的時域全顯式數(shù)值計算方法。

對本文方法進(jìn)行了必要的驗證。將該方法與透射人工邊界方法相結(jié)合，應(yīng)用于飽和兩相介質(zhì)近場波動問題，進(jìn)行了飽和土場地地震響應(yīng)的計算研究?；诒疚姆椒ㄖ袝r域遞推計算格式的傳遞矩陣，進(jìn)行了該方法穩(wěn)定性特性的研究，給出了相應(yīng)的穩(wěn)定性條件。研究結(jié)果表明：

(1)本文方法中應(yīng)用的矩陣對角化方法具有合理性。應(yīng)用本文方法得到的數(shù)值解與相應(yīng)的解析解符合良好，驗證了本文方法的正確性。

(2)基于本文方法得到的近場波動問題的計算結(jié)果符合彈性波動理論的基本規(guī)律，表明了本文方法對于典型的飽和兩相介質(zhì)近場波動問題時域計算求解的適用性。

(3)本文建立的飽和兩相介質(zhì)近場波動問題的時域數(shù)值計算方法具有全顯式算法的基本特征。方法中對動力響應(yīng)的全部分量(固相位移和速度、孔隙流體壓力)均采用顯式計算模式，即按照計算時間步數(shù)遞推和迭代的模式進(jìn)行求解，避免了求解耦聯(lián)的動力方程組，因此，該方法具有較高的計算效率，是進(jìn)行飽和兩相介質(zhì)近場波動問題時域計算求解的一種有效的算法。

本文僅是基于固相的線彈性本構(gòu)關(guān)系，開展了飽和兩相介質(zhì)彈性波動問題的研究。本文建立的時域數(shù)值計算方法對彈塑性波動問題也具有適用性。該算法與適當(dāng)?shù)膹椝苄詣恿Ρ緲?gòu)模型相結(jié)合，可以應(yīng)用于飽和兩相介質(zhì)彈塑性波動問題的研究。由于彈塑性動力本構(gòu)模型的應(yīng)力-應(yīng)變關(guān)系的表達(dá)式通常為增量形式，與此相適應(yīng)，本文算法中動力響應(yīng)的表達(dá)式也應(yīng)采用增量形式。目前相關(guān)的研究工作正在進(jìn)行。