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

(1.山東建筑大學(xué)科研處，山東 濟(jì)南 250101；2.山東建筑大學(xué)理學(xué)院，山東 濟(jì)南 250101；3.商河縣住房和城鄉(xiāng)建設(shè)服務(wù)中心，山東 濟(jì)南 251699)

0 引言

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

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

1 重心Lagrange 插值公式

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

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

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

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

根據(jù)式(4)，u(m)(x) 在節(jié)點(diǎn)xi(i ＝1，2，…，n)處的m階導(dǎo)數(shù)由式(5)表示為

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

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

2 BIOT 固結(jié)問題的重心插值配點(diǎn)法

2.1 一維線性BIOT 固結(jié)方程

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

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

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

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

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

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

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

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

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

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

2.2 BIOT 方程重心插值配點(diǎn)法計算公式

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 數(shù)值算例

MATLAB 具有強(qiáng)大的矩陣計算優(yōu)勢，算例程序均用MATLAB 語言編寫。 計算節(jié)點(diǎn)均采用區(qū)間[0，1] 上的Chebyshev 節(jié)點(diǎn)[20]，由式(29)表示為

通過坐標(biāo)變換x＝(a＋b)/2 ＋t(b-a)/2 將式(29)的節(jié)點(diǎn)變換為任意區(qū)間[a，b] 上的計算節(jié)點(diǎn)。

一維固結(jié)問題示意圖如圖1 所示。 計算參數(shù)取值為滲透系數(shù)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 一維固結(jié)問題示意圖

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

表1 不同節(jié)點(diǎn)數(shù)土體位移和孔隙水壓力最大相對誤差表

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

圖2 重心插值配點(diǎn)法計算的孔隙水壓力節(jié)點(diǎn)誤差分布圖

圖3 重心插值配點(diǎn)法計算的土體位移節(jié)點(diǎn)誤差分布圖

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

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

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

4 結(jié)論

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

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

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