重心插值配點法求解BIOT 固結問題

趙曉偉王兆清李廣惠商麗華

(1.山東建筑大學科研處，山東 濟南 250101；2.山東建筑大學理學院，山東 濟南 250101；3.商河縣住房和城鄉建設服務中心，山東 濟南 251699)

0 引言

土體固結是一個相當復雜的過程，與土體類別、邊界條件、排水條件和承載方式等有關。 BIOT[1]從彈性理論出發，研究了變形與孔隙壓力的相互作用，確保土中應力和應變滿足相容條件，給出描述比奧(BIOT)土體固結問題的偏微分方程組模型。 該模型已廣泛應用于地質力學、水文地質學、石油工程[2-3]等領域。 但是使用解析法求解BIOT 固結微分方程組很困難，因此大多采用數值法求解。 目前，常用的數值方法有限差分法、有限元法和邊界元法等。 趙維炳等[4]對BIOT 固結問題的研究采用的是中心差分格式。 但差分法對網格的規則性有要求，因此在固體力學領域受到了一定的限制。CHRISTIAN 等[5]結合有限元和有限差分法求解了BIOT 固結方程。 LEWIS 等[2]采用有限元方法研究了BIOT 理論模型，將其方程離散，再求解，并分析了地表沉降。 GASPAR 等[6-8]和BOOKER 等[9]利用有限差分方法離散了BIOT 的數學模型，并解決了多孔、分層以及二維準靜態的固結問題。 WANG等[10-11]提出了一種基于Heaviside 階梯函數的局部Petrov-Galerkin 方法以及徑向基函數的無網格方法，求解BIOT 固結模型。 有限差分法和有限元法屬于低階的數值分析方法，想要計算出更高精度的結果，有限差分法差分步長要更小，而有限元法計算網格劃分應更加稠密，這就降低了分析問題的效率。邊界元方法雖可降低問題維數，簡化問題分析，但其構造依賴于問題的基本解。 基于徑向基函數的數值方法的原理是將徑向基函數引入到未知函數的近似過程，這種方法雖然不用劃分網格，但也有其局限性，計算精度過分依賴所選區的形狀參數。

重心Lagrange 插值配點法是近年來興起的一種無網格方法，是基于微分方程強形式的配點型方法。未知函數的近似函數采用離散節點上的重心型插值。 該方法已成功應用于求解不同類型的微分方程問題，其優點是數值計算穩定性好[12]、格式較為簡單[13-14]、精度非常高[15-16]。 文章采用雙變量重心Lagrange 插值公式近似未知函數，提出了求解BIOT固結問題的重心插值配點法。 對于與時間相關的初邊值問題，傳統配點型方法在空間域上用配點格式，在時間域上用差分格式，而文章在空間域和時間域上均采用重心Lagrange 插值配點計算，實現了微分方程初邊值問題在時空域上的統一格式計算。

1 重心Lagrange 插值公式

u(x) 為定義在區間[a，b] 上的函數，其在節點a＝x1＜x2＜…＜xn＝b上的函數值為uj ＝u(xj)，j ＝1，2，…，n。u(x) 的重心Lagrange 插值由式(1)[17]表示為

重心Lagrange 插值基函數由式(2)表示為

則函數u(x) 的重心Lagrange 插值由式(3)表示為

由式(3)可知，函數u(x) 的m階導數可由式(4)表示為

根據式(4)，u(m)(x) 在節點xi(i ＝1，2，…，n)處的m階導數由式(5)表示為

將式(5)改寫成矩陣形式，由式(6)表示為

一階微分矩陣可以通過對式(2)求導直接得到，而高階微分矩陣可以通過式(7)[18]的遞推公式得到。

2 BIOT 固結問題的重心插值配點法

2.1 一維線性BIOT 固結方程

在一維固結情況下，經典BIOT 固結模型關于未知的流體壓力p(x，t) 和土體骨架位移u(x，t) 的控制偏微分方程可由式(8)表示為

式中λ、μ為土體的Lame 系數；ne為孔隙率；β為流體的壓縮系數；ke為土的滲透率；q(x，t) 為土體受力，N。

自由排水面邊界條件由式(9)表示為

剛性的不可滲透面邊界條件由式(10)表示為

初值條件由式(11)表示為

表示固結過程剛開始時水分變化為零。

將控制方程和邊界條件無量綱化，由式(12)表示為

式中f(x，t)由q(x，t) 無量綱化得到。

邊界條件由式(13)和(14)表示為

初始條件由式(15)表示為

2.2 BIOT 方程重心插值配點法計算公式

在空間域x方向和時間域t方向分別布置m和n個節點xi(i ＝1，2，…，m) 、tj(j ＝1，2，…，n) ，構成張量積型計算節點(xi，tj)(i ＝1，2，…，m，且j ＝1，2，…，n)，未知函數u(x，t) 和p(x，t) 在這些節點上的函數值為uij ＝u(xi，tj) 、pij ＝p(xi，tj) 。 其中，i ＝1，2，…，m；j ＝1，2，…，n。 利用重心Lagrange 插值公式，未知函數u(x，t) 和p(x，t) 近似由式(16)和(17)表示為

式中Li(x)、Mj(t) 分別為重心拉格朗日插值在節點xi(i ＝1，2，…，m)、tj(j ＝1，2，…，n) 上的插值基函數。

將式(16)和(17)的近似函數代入式(12)，可以得到

令式(18)在所有的計算節點(xi，tj) (i ＝1，2，…，m；j ＝1，2，…，n)上成立，得到式(19)為

注意到插值基函數性質Li(xr)＝δri、Mj(xl)＝δlj，并利用微分矩陣的記號，式(18)方程組的矩陣形式可由式(20)表示為

式中C′、D′分別為函數在節點xi(i ＝1，2，…，m) 、tj(j ＝1，2，…，n) 上重心拉格朗日插值的一階微分矩陣；C″為函數在節點xi(i ＝1，2，…，m) 上的重心拉格朗日插值的二階微分矩陣；Im、In分別為m、n階單位矩陣；符號?表示矩陣的Kroneckor 積[19]；U、P為未知函數在計算節點處函數值構成的列向量，分別由式(21)和(22)表示為

裸模這個職業由來已久，陳小北倒沒想到記者會問這個問題，他把目光投向葉曉曉，他以為她會簡單地回答自己的身體漂亮之類的話。

令L1=-(C″?In)、L2=C′?In、L3=D′?C′、L4= -ke(C″?In)＋a(Im?D′)，則式(20)可以由式(23)表示為

可進一步由式(24)表示為

2.3 邊界條件的離散和施加方法

求解式(24)的代數方程組，需要施加適當的邊界條件。 考慮邊界條件的離散問題＝-1、p＝0(x＝0，t∈(0，T))，其離散形式由式(25)表示為

式中C′1k為函數在節點xi(i ＝1，2，…，m) 上重心拉格朗日插值的一階微分矩陣的元素。

式(14)的邊界條件的離散形式由式(26)表示為

式(15)的初始條件的離散形式由式(27)表示為

式中D′1k為函數在節點tj(j ＝1，2，…，n) 上重心拉格朗日插值的一階微分矩陣元素。

把式(25)～(27)的3 個方程附加到式(24)的方程組后，可以得到新的方程組，由式(28)表示為

3 數值算例

MATLAB 具有強大的矩陣計算優勢，算例程序均用MATLAB 語言編寫。 計算節點均采用區間[0，1] 上的Chebyshev 節點[20]，由式(29)表示為

通過坐標變換x＝(a＋b)/2 ＋t(b-a)/2 將式(29)的節點變換為任意區間[a，b] 上的計算節點。

一維固結問題示意圖如圖1 所示。 計算參數取值為滲透系數ke＝1、a＝neβ(λ＋2μ)＝0、f(x，t)＝0。該算例的解析解p(x，t)＝e-0.25tcos(10π/3)sin(0.5x)和u(x，t)＝-2e-0.25tcos(10π/3)cos(0.5x)。

圖1 一維固結問題示意圖

采用重心插值配點法求解該問題，在空間和時間區域采用相同數量節點計算。 不同數量計算節點的計算結果見表1，在厚度方向和時間軸上分別取10 個節點，計算得到的土體位移和孔隙水壓力的數值解與解析解的最大相對誤差的最小值分別為1.004 7×10-14和1.215 7×10-13，計算結果接近機器精度。

表1 不同節點數土體位移和孔隙水壓力最大相對誤差表

取13 個節點時，孔隙水壓力p(x，t)數值解與精確解的相對誤差pc，e如圖2 所示；土體位移u(x，t)數值解與精確解相對誤差uc，e如圖3 所示。 在土體表面、中間和最底層處，時間從0 到1 的固結過程如圖4所示。 在表面、中間和最底處，時間從0 到20 的固結過程如圖5 所示。 由圖5 可以看出，在t＝18 時，固結過程進入穩定狀態。

圖2 重心插值配點法計算的孔隙水壓力節點誤差分布圖

圖3 重心插值配點法計算的土體位移節點誤差分布圖

圖4 土體固結過程圖(t 從0 到1 )

圖5 土體固結過程圖(t 從0 到20 )

該算例在選取節點數較少時，其位移和孔隙水壓力的數值解與精確解的相對誤差可達10-13～10-14，已經接近MATLAB 的機器精度，數值計算公式的正確性和計算方法的有效性由此可以得到驗證，表明了該方法具有極高的計算效率。 當節點數量持續增加時，數值結果變化較小，表明該方法具有極好的數值穩定性。

4 結論

文章采用雙變量重心Lagrange 插值公式近似未知函數，提出了求解BIOT 固結問題的重心插值配點法，并針對土體結構的BIOT 固結方程， 并通過算例驗證，得到的主要結論如下:

(1) 采用重心插值配點法離散控制方程，不需積分計算，也不用劃分網格，極大地提高了計算效率。 離散點的分布采用切比雪夫節點，可使計算結果具有良好的數值穩定性。 數值算例表明，所用重心插值配點法沒有出現runge 現象，即使增加計算節點，計算結果依然具有穩定的、高精度的收斂特性，其計算誤差非常接近機器精度。

(2) 所用重心插值配點法的計算公式用矩陣-向量形式表達，可使計算程序的編寫過程變得簡單，易于理解。 邊界條件的施加采用附加法，將離散控制方程和邊界條件約束的代數方程組合起來，任意的邊界條件都可以用該方法實現。