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

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

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

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

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

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

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

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

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

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

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

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

向量 m的表達式為：

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

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

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

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

2 方法的建立

2.1 波動方程系數矩陣對角化

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

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

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

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

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

2.2 動力響應的時域顯式逐步計算列式

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

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

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

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

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

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

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

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

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

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

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

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

在迭代計算的每個時間分步上，動力響應的求解須按照如下步驟進行：1)分別應用式(23)和式(25)由前一個時間分步的動力響應計算當前分步的固相位移uk+1和孔隙流體壓力pk+1；2)應用式(29)由前一時間分步的動力響應以及步驟1)得到的uk+1和 pk+1計算當前分步的固相速度u˙k+1。按照計算時間步數對上述求解步驟進行循環，即可完成飽和兩相介質動力響應時程的計算。基于本文算法的求解模式可知，飽和兩相介質波動問題相關的全部動力響應，即固相位移 u 和固相速度 u˙，以及孔隙流體壓力 p都是采用顯式的計算模式，即迭代和遞推的模式進行時域求解，因此本文建立的算法是一種全顯式的時域數值求解方法。

3 方法的驗證

3.1 矩陣對角化合理性驗證

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

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

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

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

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

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

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

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

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

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

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

3.2 方法的解析解驗證

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

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

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

圖7 飽和土柱動力響應計算模型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 飽和土柱材料參數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 可知，應用本文算法得到的飽和土柱動力響應的數值解結果與相應的解析解結果符合良好，表明本文算法能夠對飽和兩相介質的動力響應進行較為準確的計算，驗證了本文算法的正確性。

4 方法的應用

本節將本文建立的時域顯式數值計算方法應用于飽和兩相介質的近場波動問題之中，進行飽和土場地地震響應的計算研究。本項研究的意義表現在兩個方面：第一，場地的地震響應特性是土動力學和巖土地震工程領域中的重要研究課題之一[25]；第二，場地地震響應的計算是地震動輸入條件下場地土體與結構動力相互作用問題研究的基礎和重要步驟。

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

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

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

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

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

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

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

圖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 方法穩定性研究

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

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

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

令

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

式中：

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

式中：

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

以時間步長表示的穩定性條件為：

6 結論

本文基于以固相位移u 和孔隙流體壓力p 作為基本未知量的u-p 形式的飽和兩相介質彈性波動方程，開展了飽和兩相介質近場波動問題時域顯式數值計算方法的研究。通過對波動方程中的質量矩陣和孔隙流體壓縮矩陣進行對角化處理，消除了方程中的動力耦聯，實現了波動方程的解耦。基于解耦的有限元波動方程組，分別應用中心差分法和Newmark 常平均加速度法求解固相位移和速度，基于向后差分法求解孔隙流體壓力，推導得到了飽和兩相介質動力響應的時域顯式逐步積分計算的計算列式，建立了飽和兩相介質近場波動問題的一種新的時域全顯式數值計算方法。

對本文方法進行了必要的驗證。將該方法與透射人工邊界方法相結合，應用于飽和兩相介質近場波動問題，進行了飽和土場地地震響應的計算研究。基于本文方法中時域遞推計算格式的傳遞矩陣，進行了該方法穩定性特性的研究，給出了相應的穩定性條件。研究結果表明：

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

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

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

本文僅是基于固相的線彈性本構關系，開展了飽和兩相介質彈性波動問題的研究。本文建立的時域數值計算方法對彈塑性波動問題也具有適用性。該算法與適當的彈塑性動力本構模型相結合，可以應用于飽和兩相介質彈塑性波動問題的研究。由于彈塑性動力本構模型的應力-應變關系的表達式通常為增量形式，與此相適應，本文算法中動力響應的表達式也應采用增量形式。目前相關的研究工作正在進行。