基于分形理論的飽和土孔隙水壓力計算模型

2楊忠翰田東弘

(1.東北林業(yè)大學(xué) 土木工程學(xué)院， 哈爾濱 150040；2.濟寧市勘測院, 山東 濟寧 272000 )

飽和土作為一種兩相介質(zhì)，是由土顆粒及孔隙中的流體混合而成的。因此在研究飽和土體的應(yīng)力-應(yīng)變關(guān)系時，必須要考慮土顆粒與流體間的相互作用關(guān)系。在受到外界應(yīng)力作用時，土顆粒骨架的壓縮會使飽和土體中孔隙水壓力升高，導(dǎo)致土體的有效應(yīng)力降低，并使得土體的抗剪強度下降，從而造成土體的液化。對于孔隙水壓力的計算，土顆粒間孔隙中的流體流動及流量變化對孔隙水壓力的影響往往被忽略，在計算孔隙水壓力時，計算值往往小于實際值，導(dǎo)致人們對于土體的有效應(yīng)力及抗剪強度計算存在誤差。本文基于孔隙數(shù)目-尺寸分形模型[1-2]，推導(dǎo)出孔隙度與分形維數(shù)之間的數(shù)學(xué)關(guān)系，并進一步基于飽和土中孔隙水流動方程及水壓力方程，推導(dǎo)出孔隙間流量變化導(dǎo)致的孔隙水壓力增量與分形維數(shù)及流量變化的定量關(guān)系，并進一步導(dǎo)出飽和土體中任一點孔隙水壓力計算公式。

1 巖土材料的孔隙度與分形維數(shù)的演化關(guān)系

孔隙-固體分形模型形成于20世紀90年代，早在1989年，Neimark[3]發(fā)展了自相似多尺度逾滲系統(tǒng)，用于表征孔-固相界面為分形的無序離散的多孔介質(zhì)；1994年，綜合分形孔隙和固體的數(shù)量-尺寸分布，Perrier[4]獨立地提出了一個土壤結(jié)構(gòu)的多尺度模型；基于Neimark給出的定義，Perrier等[5]進一步系統(tǒng)地建立了孔隙-固體分形模型。該模型可用于表征具有表面分形、分形孔隙尺寸分布、分形固體質(zhì)量及分形孔隙質(zhì)量特性的多孔介質(zhì)。因為巖土類材料從微觀角度講，是一種由大小各不相同的土顆粒堆積而成的一種多孔介質(zhì)，此類材料的構(gòu)造過程與分形十分相似，因此，可以基于分形模型對巖土類材料的孔隙度進行定量分析。

1.1 分形維數(shù)的計算

基于孔隙-固體分形模型給出的定義，在一個自相似系統(tǒng)中，分形物體的累計測量數(shù)目與單元尺寸之間通常存在冪律關(guān)系[6]，即

N(ri)=k(ri)-D,i=0,1,2,...,。

(1)

則分形維數(shù)D為

(2)

式中：N(ri)是長度為ri的單元個數(shù)；k為單位長度的初始元個數(shù)；i為迭代次數(shù)；D為分形維數(shù)。式(1)為基于分形理論的孔隙和固體的累計數(shù)目-尺寸分布的普適關(guān)系，與自相似系統(tǒng)中初始元維數(shù)及構(gòu)造方法無關(guān)。

1.2 孔隙度的計算

1992年，Tyler和Wheatcraft曾指出土體材料中可能同時存在孔隙數(shù)目-尺寸和固體數(shù)目-尺寸分形冪律關(guān)系[7]，但當(dāng)時還沒有解釋此現(xiàn)象的理論。而Perrier等[5]所提出的孔隙-固體分形模型則提供了一種明確的理論模型。鄭瑛等[8]基于此模型對多孔材料的孔隙結(jié)構(gòu)進行了分形描述。本文基于孔隙-固體分形模型的構(gòu)造過程，對土體材料孔隙度與分形維數(shù)之間的關(guān)系進行如下推導(dǎo)。

Vs=(n3-m)i(R/ni)3。

(3)

則孔隙體積可表示為

Vp=R3-(n3-m)i(R/ni)3=

(4)

則孔隙度可表示為

(5)

由式(2)可得出立方體分形維數(shù)為

(6)

式中：n的直觀含義為初始元樣本內(nèi)某一邊長方向上的樣本邊長與最高階次孔隙平均直徑之比，即樣本的測量尺度(初始元等份數(shù))；m表示最高階次孔隙的個數(shù)[3]。

由此可見，分形維數(shù)D可以反映樣本內(nèi)孔隙的大小及其分布空間的關(guān)系。

由式(6)可得，m=n3-nD，代入式(5)可得孔隙度與分形維數(shù)的數(shù)學(xué)關(guān)系式，即

φ=1-(nD-3)i。

(7)

該估計是先驗過程的期望分布F0(x)和經(jīng)驗分布估計Fn(x|x1,x2,···,xn)的加權(quán)平均.由于經(jīng)驗分布函數(shù)是階梯函數(shù),為得到一個光滑的分布估計,用核估計代替經(jīng)驗分布函數(shù)Fn(x|x1,x2,···,xn),則總體的密度估計為

(8)

由式(8)可以看出，在測量尺度一定的情況下，通過分形維數(shù)即可求得土體多孔介質(zhì)的孔隙度?；诖斯剑蛇M一步利用分形維數(shù)對土體孔隙水壓力進行定量描述。

2 飽和土孔隙水壓力計算模型

2.1 測量球概念的提出

對于完整的巖土材料而言，孔隙水存在于其內(nèi)部的土顆粒間的孔隙當(dāng)中。想要通過模型反映土體內(nèi)部孔隙水壓力的大小及分布情況，準確地表征土體內(nèi)部孔隙的分布情況是十分重要的。因此，在考慮孔隙水壓力時，Okada等[9]提出了測量球的概念。測量球是一個虛擬球體空間，其球心與土顆粒球心重合，半徑為土顆粒半徑2倍。如圖1所示，測量球是一個由若干土顆粒及顆粒間流體組成的兩相體系，其優(yōu)點在于可以將飽和土樣的孔隙度與孔隙水壓力等物理及力學(xué)參數(shù)作離散處理。根據(jù)Okada等的研究，測量球半徑為對應(yīng)的顆粒半徑2倍時，空間內(nèi)部包含土體顆粒數(shù)目適中，能夠獲得準確的體積應(yīng)變。該模型假定固體顆粒為剛體，而孔隙中流體的體積可以發(fā)生改變，且流體的體積變化可以產(chǎn)生對應(yīng)的孔隙水壓力。

圖1 Okada的測量球示意圖Fig.1 Okada’s measuring sphere

2.2 孔隙水壓力增量計算

飽和土體為固液兩相介質(zhì)，土顆粒間的孔隙全部充滿流體，因此對于飽和土中任意位置測量球，其孔隙度φ可以定義為

(9)

式中：φ表示測量球孔隙度；Vv表示測量球孔隙體積；Vw表示流體體積；V表示測量球體積。將式(7)代入式(9)得

Vw=φV=[1-(nD-3)i]V。

(10)

圖2表示相鄰兩測量球的位置關(guān)系，假設(shè)測量球內(nèi)的土顆粒與孔隙為均勻分布，故兩測量球的過流面積可表示為兩測量球的接觸面積與平均孔隙度乘積，由Darcy定律可知，單位時間內(nèi)流過兩測量球交界面的滲透流量為

(11)

圖2 相鄰兩測量球位置關(guān)系Fig.2 Position relation between two adjacent measuring spheres

則A，B兩測量球之間Δt時間內(nèi)的流量可表示為

(12)

式中：pA，pB表示兩測量球的孔隙水壓力；ρ表示20 ℃水的密度；g表示重力加速度。由此可知，孔隙度可以確定相鄰兩測量球的過流面積，進而確定孔隙水流的交換速度。

本文基于相鄰兩測量球中心距離與測量球半徑的數(shù)值關(guān)系，分別計算出不同情況下兩測量球的接觸面積，即過流面積。相鄰兩測量球的位置關(guān)系可分為相交或相離2類。其中相交情況又可分為兩測量球部分相交或半徑較小測量球處于半徑較大測量球之內(nèi)。當(dāng)兩測量球相交時，接觸面積可通過兩測量球A,B的相對位置進行計算。不妨假設(shè)測量球A與測量球B的半徑關(guān)系為RA

其中Rt=0.5(RA+RB+LAB)。當(dāng)兩測量球相離，即LAB≥RA+RB時，S=0；當(dāng)半徑較小測量球處于半徑較大測量球之內(nèi)，即LAB≤RB時，S=πRB2。

在Δt時間內(nèi)，一個測量球與其周圍相鄰測量球發(fā)生孔隙水流量交換，假設(shè)在流量交換過程中，測量球孔隙度不發(fā)生改變，則交換過程中，測量球的流量總變化可表示為∑ΔQ。根據(jù)孔隙水壓力方程可得流量交換導(dǎo)致的孔隙水壓力增量，即

(14)

式中Ew為孔隙水壓縮模量。將式(10)代入式(14)可得

(15)

式(15)可視為飽和狀態(tài)下孔隙水壓力增量與土體分形維數(shù)之間的函數(shù)關(guān)系。實際上在飽和土體內(nèi)部發(fā)生孔隙水流量交換時，土體孔隙度及滲透系數(shù)會發(fā)生微小的變化。此處假定孔隙度不變且滲透系數(shù)為常數(shù)是作簡化處理。

2.3 測量球內(nèi)土顆粒所受水壓力計算

土體在飽和狀態(tài)下，其內(nèi)部完全由土顆粒及孔隙間流體組成。其中對于完全浸沒在流體中的土顆粒而言，所受浮力遠小于孔隙水壓力，故浮力可忽略不計。假設(shè)測量球內(nèi)各孔隙均為相互連通狀態(tài)，即測量球內(nèi)部孔隙水壓力處處相同。如圖3所示，只考慮單一測量球內(nèi)部孔隙水壓力的情況下，完全浸沒在測量球內(nèi)部流體中的土顆粒所受流體壓力是平衡的；而邊界處部分浸沒的土顆粒所受流體壓力為非平衡的。

圖3 測量球內(nèi)不同位置孔隙水壓力示意圖Fig.3 Schematic diagram of pore water pressure at different locations in the measuring sphere

因此，對于測量球內(nèi)部任意一點的孔隙水壓力，只需考慮邊界處位置即可。對于邊界區(qū)域的寬度，應(yīng)滿足測量球內(nèi)最大顆粒位于邊界區(qū)域時，該顆粒部分浸沒于孔隙流體中而非完全浸沒。因此可根據(jù)最大顆粒半徑限定邊界區(qū)域?qū)挾?，?

式中：x,y,z分別表示測量球球心的三維坐標(biāo)；xi,yi,zi分別表示測量球內(nèi)任一土顆粒中心的三維坐標(biāo)。

確定測量球邊界條件后，將孔隙水壓力在其作用面上積分，即可得到測量球不平衡區(qū)域內(nèi)任一點孔隙水壓力，即

(17)

式中：ri為邊界土顆粒i與測量球接觸區(qū)域所構(gòu)成圓的半徑。式(18)即孔隙水在飽和土孔隙間發(fā)生流量交換時，假設(shè)土體孔隙度不發(fā)生改變的前提下，推導(dǎo)的基于分形理論的飽和土中任意一點孔隙水壓力計算模型。

3 數(shù)值算例

選取一個飽和的均質(zhì)粉質(zhì)黏土邊坡作為算例，以驗證本文所建立的孔隙水壓力模型實用性及有效性。算例的邊坡幾何模型即網(wǎng)格如圖4所示。模型全長為1 610 m，高505 m，其中邊坡長680 m，邊坡高度306.5 m。模型所選飽和粉質(zhì)黏土的主要力學(xué)參數(shù)詳見表1。基于本文所建立的孔隙水壓力模型，算例中孔隙水壓力設(shè)定為與粉質(zhì)黏土分形維數(shù)相關(guān)的函數(shù)。算例中邊坡內(nèi)部孔隙水流速變化很小，因此可假設(shè)單位時間內(nèi)土體中孔隙水流量變化為一常數(shù)ΔQ，則邊坡中任一高度h處孔隙水壓力可表示為

(19)

式中粉質(zhì)黏土相關(guān)參數(shù)n和D可從土壤的水分特征曲線數(shù)據(jù)的擬合結(jié)果[10]中得出，代入后即可在粉質(zhì)黏土邊坡不同高度施加相應(yīng)的孔隙水壓力。

圖4 算例幾何模型Fig.4 Geometric model of the computation example

彈性模量/MPa泊松比內(nèi)摩擦角/(°)黏聚力/kPa重度/(kN·m-3)1 5000.35343323

圖5為算例中邊坡飽水狀態(tài)下的塑性區(qū)位移分布，其位移分布規(guī)律基本與文獻[11]中相似條件下的位移分布規(guī)律一致。邊坡失穩(wěn)通常發(fā)生在降雨過程中，隨土中含水率升高，土體逐漸趨于飽和，孔隙水壓力上升致使土體有效應(yīng)力降低，抗剪強度下降。邊坡中飽水的土體區(qū)域持續(xù)產(chǎn)生塑性變形，區(qū)域內(nèi)不平衡力范圍逐漸擴大，最終會導(dǎo)致邊坡失穩(wěn)產(chǎn)生滑坡。

圖5 塑性區(qū)位移分布Fig.5 Contours of plastic zone displacement

本文中孔隙水壓力計算模型考慮飽和土中孔隙水滲流時流量交換產(chǎn)生的孔隙水壓力增量，與以往基于傳統(tǒng)孔隙水壓力計算公式所做的邊坡穩(wěn)定分析相比，更能反映邊坡產(chǎn)生滑動破壞的真實情況?；诜中卫碚摰玫降目紫端畨毫τ嬎愎?，可通過不同類型土體對應(yīng)的分形維數(shù)研究不同工況中各種土體內(nèi)部孔隙水壓力的數(shù)值及分布情況，對實際工程具有重要的參考價值。

4 結(jié) 語

基于孔隙數(shù)目-尺寸分形模型建立了土體材料孔隙度與分形維數(shù)之間的關(guān)系，并且進一步推導(dǎo)了飽和土體孔隙水壓力狀態(tài)方程，建立了飽和土體中任一位置土顆粒所受孔隙水壓力與分形維數(shù)及孔隙間流量變化之間的關(guān)系。并將所建立的飽和土孔隙水壓力函數(shù)代入ANSYS命令流中，對飽和狀態(tài)下的邊坡進行有限元穩(wěn)定性分析，以驗證建立的飽和土孔隙水壓力計算模型的準確性及實用性。本文所得公式可用于飽和土體的有效應(yīng)力及抗剪強度計算修正，并可應(yīng)用于飽和土體宏觀-微觀的多尺度液相-固相耦合滲流分析。