張 昭，程靖軒，劉奉銀，齊吉琳，許增光

（1.西安理工大學(xué) 巖土工程研究所，陜西 西安 710048；2.北京建筑大學(xué) 土木與交通工程學(xué)院，北京 100044；3.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國家重點實驗室，陜西 西安 710048）

1 研究背景

土壩常因壩前水位驟降而造成壩體滑坡，這種現(xiàn)象與非飽和滲流對土坡變形及其穩(wěn)定性的影響密切相關(guān)［1-2］。土壩在水位驟降時滲流前鋒邊緣的后退主要受土顆粒間的毛細黏聚作用影響［3］，因此，在采用本構(gòu)模型進行土壩變形分析及其穩(wěn)定性評價時，如何考慮這種毛細黏聚作用的細觀機理值得深入研究。當壩前水位驟降時，在構(gòu)成土壩的濕顆粒材料中，水分的遷移使其水力與力學(xué)行為發(fā)生變化，兩個表面濕潤的顆粒相互靠近時水分會在其接觸點及附近區(qū)域形成液橋，使顆粒發(fā)生連接［4］，這種現(xiàn)象取決于毛細黏聚作用［5-6］，能夠從細觀角度反映水分含量對濕顆粒材料變形和強度的貢獻［7-8］。若要定量描述濕顆粒材料的毛細黏聚作用機制，則可將其簡化為一對不等徑球體濕顆粒，采用液橋描述濕顆粒間吸附的水分形態(tài)，通過分析液橋的幾何形狀與受力形態(tài)來構(gòu)建液橋斷裂距離及其毛細力關(guān)于其體積和顆粒間距的關(guān)系式［9-12］。由圖1可知，在濕顆粒材料因水位驟降從完全飽和狀態(tài)開始減濕的過程中，水分會經(jīng)歷毛細狀態(tài)（顆粒被水分完全浸沒）、索帶狀態(tài)（形如索帶狀液橋的水分至少同時與3 個顆粒相連）和鐘擺狀態(tài)（形如鐘擺狀液橋的水分僅與兩個顆粒相連）［5，13］。本文僅以鐘擺狀液橋為研究對象，并認為其表面關(guān)于水平軸對稱，處于擬靜力平衡狀態(tài)，不考慮黏滯作用和重力作用［14］對其力學(xué)特征的影響。

圖1 壩前水位驟降時濕顆粒與水的相互作用原理

液橋的力學(xué)特征取決于其表面幾何形狀［9，11-12，15］，已有學(xué)者通過濕顆粒間液橋拉伸試驗研究［10，16-20］表明Young-Laplace 方程能準確描述液橋表面幾何形狀，并可認為其表面平均曲率為常數(shù)。求解Young-Laplace 方程可得到液橋的幾何特征、斷裂距離和毛細力［9，11-12，15］。

學(xué)者們主要采用以下4種方法研究濕顆粒間液橋的斷裂距離及其毛細力：①液橋拉伸試驗［10，16-17，19，21］：雖可直接測定液橋斷裂距離及其毛細力，卻需在試驗過程中避免重力和顆粒表面粗糙度的影響；②Young-Laplace 方程的數(shù)值解［9-12］：只要對這種非線性方程進行合理地離散化，則求取的數(shù)值解在允許誤差范圍內(nèi)較為準確；③曲線擬合公式［10-11］：依據(jù)液橋拉伸試驗與Young-Laplace 方程數(shù)值解的數(shù)據(jù)組構(gòu)建，雖能在較大的液橋體積與顆粒間距范圍內(nèi)實現(xiàn)快速、準確計算，但需標定擬合參數(shù)；④半解析-半數(shù)值方法［12，22-28］：對液橋表面形狀進行簡化假定，故具有明確的物理意義，在液橋體積與顆粒間距的一定范圍內(nèi)預(yù)測較為準確。

此外，依據(jù)上述4 種方法的研究成果亦可知：顆粒粒徑對液橋的幾何特征與受力狀態(tài)影響顯著。一方面，當采用大、小顆粒半徑之比（下文簡稱顆粒半徑比）描述這種影響時，液橋的毛細力會隨顆粒半徑比增大呈遞增趨勢［10，16］。另一方面，在曲線擬合與半解析-半數(shù)值方法中，常采用顆粒等效半徑將等徑濕顆粒間液橋的斷裂距離及其毛細力公式推廣至不等徑濕顆粒的情況［10-11］中，顆粒等效半徑Re可用大、小顆粒半徑（R1和R2）表示為：

式中：當濕顆粒半徑相同時，R1=R2=R（mm），則Re=R。值得注意：Willett 等［10］依據(jù)Young-Laplace方程數(shù)值解分析了采用顆粒等效半徑描述不等徑濕顆粒間液橋斷裂及其毛細力變化特征的準確性，結(jié)果發(fā)現(xiàn)：采用顆粒等效半徑預(yù)測的小體積液橋毛細力誤差小于2%，而對大體積液橋毛細力的預(yù)測誤差高達20%；此外，Willett 等［10］發(fā)現(xiàn)采用顆粒等效半徑預(yù)測的液橋斷裂距離相對偏大，進而對相應(yīng)的Young-Laplace 方程數(shù)值解進行曲線擬合構(gòu)建了液橋斷裂距離關(guān)于其體積、固-液接觸角和顆粒半徑比的擬合公式。類似地，Lian 等［11］也對Young-Laplace 方程數(shù)值解進行曲線擬合，雖得到了適用于小體積液橋斷裂距離及其毛細力的計算公式，卻因未考慮顆粒半徑比的影響而不適用于不等徑濕顆粒間的大體積液橋。

不僅如此，半解析-半數(shù)值方法大都假定液橋表面在二維平面內(nèi)形如圓?。磮A弧假定）［22-24，26-28］。該假定使用簡便，但毛細力在液橋最窄“頸部”和顆粒-液橋接觸點處不同的現(xiàn)象與Young-Laplace 方程中液橋表面的平均曲率及其毛細力均為常數(shù)的性質(zhì)不符［9，12］。為此，Kruyt 等［12］假定等徑濕顆粒間液橋表面在二維平面內(nèi)形如橢圓弧（即橢圓弧假定），不僅考慮了Young-Laplace 方程中液橋表面平均曲率為常數(shù)的性質(zhì)，而且無需引入擬合參數(shù)。然而，Kruyt 等［12］僅以等徑濕顆粒間的小體積液橋為研究對象，未考慮不等徑顆粒間的大、小體積液橋情況。

為此，以濕顆粒內(nèi)形如液橋的水分形態(tài)為研究對象，從顆粒粒徑和液橋體積這兩個影響因素出發(fā)，將解析表達式推導(dǎo)與Young-Laplace 方程數(shù)值解的數(shù)據(jù)擬合分析相結(jié)合，力圖提出不等徑濕顆粒間的大、小體積液橋斷裂距離及其毛細力的計算方法。首先，將文獻［12］中針對等徑濕顆粒間小體積液橋基于橢圓弧假定構(gòu)建的解析公式推廣至不等徑濕顆粒的范圍；其次，分析顆粒等效半徑在不同顆粒半徑比、液橋體積、顆粒間距和固-液接觸角條件下的適用性；最后，構(gòu)建考慮顆粒半徑比的斷裂距離擬合公式，并改進已有的毛細力擬合公式，使之也適用于大體積液橋，同時，將已有的與改進的擬合公式預(yù)測結(jié)果進行對比分析。本文提出的毛細力計算方法可為水位驟降時土壩的變形分析及其穩(wěn)定性評價提供一種綜合考慮細觀毛細黏聚作用機理與飽和度變化過程的壩體材料本構(gòu)計算途徑。

2 濕顆粒間液橋的幾何形狀及其Young-Laplace 方程的數(shù)值解

本節(jié)先描述不等徑濕顆粒間鐘擺狀液橋的幾何形狀及控制其表面形狀的Young-Laplace 方程，進而提出針對不等徑濕顆粒的Young-Laplace 方程數(shù)值求解方法，構(gòu)建高精度數(shù)值解的數(shù)據(jù)組，不僅可驗證橢圓弧假定，而且可用于分析不同顆粒半徑比的濕顆粒在不同顆粒間距下形成大、小體積液橋的斷裂及其毛細力變化特征。

2.1 液橋的幾何形狀從理論上分析不等徑濕顆粒間鐘擺狀液橋應(yīng)滿足的幾何條件。先將濕顆粒簡化為圖2所示一對半徑分別為R1和R2（R1≥R2，下標1 和2 分別表示大、小顆粒）的球體顆粒，不考慮顆粒的重力和浮力作用，球體顆粒被液橋不完全浸潤，故需考慮固-液接觸角θ。依據(jù)文獻［29］，θ值主要受濕顆粒的礦物成分影響，而在進行土壩的滲流與變形計算時常選取同種濕顆粒材料構(gòu)筑壩體［1，3］，因此本文認為這對不等徑濕顆粒的礦物成分相同，則液橋?qū)Υ?、小顆粒的θ值相同。關(guān)于液橋的表面形狀，當其體積較小時（VL?B≤1×10-3，VL?B為無量綱液橋體積，詳見2.3 節(jié)），可將其假定為橢圓??；當其體積較大時（VL?B＞1×10-3），可采用Young-Laplace 方程數(shù)值解描述，詳見第3-4 節(jié)。這對不等徑濕顆粒與液橋的位置關(guān)系亦如圖2所示：大、小顆粒的中心分別為O1和O2，β1和β2分別為液橋?qū)Υ?、小顆粒的充填角；P1(xc1，yc1)和P2(xc2，yc2)分別為大、小顆粒與液橋的接觸點；液橋關(guān)于x軸對稱，其表面形狀可用函數(shù)y(x)描述，實質(zhì)表示液橋在對稱軸上任意位置x處的局部半徑y(tǒng)，以其y≥0 的上半部分為研究對象，令液橋最窄“頸部”（呈圓形截面）位于x=0 處，此處y的最小值稱為液橋的最窄“頸部”半徑y(tǒng)0；大、小顆粒至x=0 處的間距分別為d1和d2。由圖2亦可知，液橋表面形狀函數(shù)y(x)在其最窄“頸部”處滿足如下條件：

圖2 不等徑濕顆粒鐘擺狀液橋的幾何參數(shù)及其受力狀態(tài)

式中：y′(x)=dy/dx。接觸點P1和P2的坐標可分別表示為：

當兩個顆粒-液橋接觸半徑（yc1和yc2）較小時（yc1/R1?1 且yc2/R2?1），式（3）可近似為：

此外，依據(jù)文獻［28］，y(x)在兩個接觸點處的斜率（y′(xc1)和y′(xc2)）與θ有關(guān)（圖2）：

式中：ξ=tanθ。需注意：本文僅討論凹形鐘擺狀液橋，因此β1和β2與θ需滿足β1+θ≤90°和β2+θ≤90°。

2.2 Young-Laplace 方程與液橋的受力狀態(tài)依據(jù)文獻［10，16-20］，液橋表面形狀可用Young-La?place 方程描述，依據(jù)該方程可認為液橋表面平均曲率為常數(shù)。本文以軸對稱液橋為研究對象，則表征其表面形狀的Young-Laplace 方程為：

式中：K[y]為平均曲率算子；y′=dy/dx；；Cm為平均曲率，m-1；ψ表示基質(zhì)吸力，kPa，即橫貫于氣-液交界面內(nèi)外的壓力差；σw為水分的表面張力（20 ℃時純水的σw=0.0728 N/m）。該非線性二階常微分方程的邊界條件如式（2）、式（3）和式（5）所示。

對式（6）所示Young-Laplace 方程求積分可得［12］：

式中：λ為毛細長度系數(shù)，mm。Λ[y]為對K[y]的第一次積分結(jié)果。式（7）描述了液橋沿其長度方向的受力平衡條件［11-12］。將式（2）和式（5）代入式（7），則Cm和λ可表示為：

關(guān)于液橋的受力狀態(tài)，選取液橋的最窄“頸部”截面進行分析，亦如圖2所示，相應(yīng)的毛細力F（mN）由作用在該截面周的表面張力σw和截面內(nèi)的基質(zhì)吸力ψ共同產(chǎn)生：

將式（6）所示ψ與σw的關(guān)系式和式（8）所示λ的公式代入式（9），即可進一步將F表示為：

由式（8）可知，y0與yc1和yc2之間存在聯(lián)系。將式（5）代入式（8）可得這三者之間的關(guān)系式：

式中：系數(shù)δ4、δ2和δ0可用y0、yc2、R1、R2和ξ表示為：

當兩個顆粒半徑相同時（R1=R2），式（12）有唯一解yc1=yc2，此時液橋也關(guān)于y軸對稱。表征式（12）所示方程存在實解的判別公式為Δ=δ22-4δ4δ0。該方程存在兩個實解（Δ＞0）、一個實解（Δ=0）或無實解（Δ＜0）均需滿足條件yc2≤yc1≤R1。當存在兩個實解時，僅液橋表面積較小（其表面能較?。┑慕庥形锢硪饬x［12］。式（12）關(guān)于大顆粒-液橋接觸半徑y(tǒng)c1存在實解時液橋最窄“頸部”半徑y(tǒng)0和小顆粒-液橋接觸半徑y(tǒng)c2的取值區(qū)域如圖3所示，實質(zhì)上給出了R1/R2= 1、2 和128 且θ=0°和20°時的上、下限。

當yc2和y0較小時，方程判別式Δ為零，且Δ可近似表示為：

式（14）可在yc2和y0處進行泰勒展開，而與顆粒半徑比無關(guān)。當ξ給定時，求解式（14）可得的最小值，即圖3中下邊界在yc2=0 處的斜率。由圖3可知，當較小時且大于1 時，下邊界基本不受影響。需注意：當較大時，與θ會顯著影響式（12）關(guān)于yc1存在實解的區(qū)域，故不同與θ條件下式（12）關(guān)于yc1存在實解區(qū)域的顯式公式有待進一步研究。

2.3 Young-Laplace 方程的數(shù)值解式（6）所示Young-Laplace 方程的數(shù)值解可依據(jù)如下方法求取。當給定兩個顆粒-液橋接觸半徑（yc1和yc2）與液橋的最窄“頸部”半徑y(tǒng)0時，平均曲率Cm和毛細力系數(shù)λ可由式（8）確定。由式（7）可知，y′=dy/dx可表示為y的函數(shù)，則液橋表面形狀函數(shù)y(x)可通過dx/dy=1/(dy/dx)求得其反函數(shù)x(y)，進而將Kruyt 等［12］計算等徑濕顆粒間液橋表面形狀x(y)的積分公式推廣至不等徑濕顆粒。該函數(shù)在x＜0 和x≥0 時的斜率（dx/dy）分別為負值和正值（圖2），故可分為x1(y)（x＜0）和x2(y)（x≥0）兩部分，依據(jù)式（7）得到x1(y)和x2(y)的積分公式：

圖3 式（12）關(guān)于yc1 存在實解時y0 和yc2 的取值區(qū)域

由于Kruyt 等［12］采用的Gauss-Kronrod 自適應(yīng)數(shù)值積分算法可解決式（15）在y0處積分出現(xiàn)的弱奇點問題，故采用該算法（令計算中的相對誤差限為10-6）計算式（15）中的積分，即可得到y(tǒng)(x)。從大、小顆粒與液橋的接觸面至液橋最窄“頸部”處（x=0）分別對應(yīng)的液橋體積（VLB1和VLB2）可通過對y(x)進行數(shù)值積分得到：

式中：VCA1和VCA2分別為大、小顆粒被液橋浸潤的球冠體積，可表示為：

液橋體積VLB（μL）及顆粒間距2d（mm）可由大、小顆粒對應(yīng)的液橋體積及顆粒間距（d1和d2）求和而得：

由此可依據(jù)顆粒等效半徑Re定義無量綱液橋體積VL?B（VL?B =VLB/Re3）和無量綱顆粒間距2d?（2d?=2d/Re）。

因此，可依據(jù)本節(jié)上述計算方法對Young-Laplace 方程求得下列條件下的不同數(shù)值解：①在充填角β低于60°的范圍內(nèi)取500 個不同數(shù)值；②液橋最窄“頸部”半徑及其與小顆粒接觸半徑之比y02在0～1 范圍內(nèi)取500 個不同數(shù)值；③顆粒半徑比R1R2=1、4/3、2、4、8、16 和128；④固-液接觸角θ=0°、20°、40°。

將上述條件下的不同數(shù)值解構(gòu)成數(shù)據(jù)組，要求在y02=0～1 的范圍內(nèi)僅考慮如圖2所示的凹形鐘擺狀液橋，而且忽略顆粒間距2d為負值時的無物理意義數(shù)值解（2d可通過y0、yc1和yc2算得）。由于該數(shù)據(jù)組存儲了每個β和y0yc2對應(yīng)的2d、VLB、F和Cm，因此，一方面可將Cm代入式（6）算得相應(yīng)的ψ，另一方面可依據(jù)Kruyt 等［12］提出的數(shù)值線性內(nèi)插方法算得VLB和2d任意取值組合對應(yīng)的F，從而為驗證、評價本文以下提出的液橋斷裂距離及其毛細力的解析與擬合公式提供數(shù)據(jù)支撐。

3 基于橢圓弧假定計算不等徑濕顆粒間小體積液橋的斷裂距離及其毛細力

所謂橢圓弧假定，即在二維平面內(nèi)采用橢圓弧描述濕顆粒間液橋的表面形狀，該假定已用于研究等徑濕顆粒間小體積液橋的斷裂及毛細力變化特征［12］。本節(jié)將橢圓弧假定推廣至不等徑濕顆粒間形成的液橋：先闡述如何采用橢圓方程描述不等徑濕顆粒間液橋的表面形狀，以討論橢圓弧假定的適用性；繼而取液橋內(nèi)毛細力為常數(shù)，構(gòu)建相應(yīng)的幾何關(guān)系式，從而確定基于橢圓弧描述液橋表面形狀所需的幾何參數(shù)；再推導(dǎo)顆粒-液橋接觸半徑、液橋斷裂距離及其毛細力關(guān)于液橋體積和顆粒間距的解析公式，通過對比分析解析公式對不同液橋斷裂距離及其毛細力的Young-Laplace 方程數(shù)值解的預(yù)測效果，旨在驗證采用顆粒等效半徑將等徑濕顆粒間小體積液橋的斷裂距離及其毛細力解析公式推廣至不等徑濕顆粒范圍的適用性。

3.1 基于橢圓方程描述不等徑濕顆粒間液橋的表面形狀橢圓弧假定可用于描述等徑濕顆粒間液橋

的表面形狀［12］。然而，對不等徑濕顆粒間的液橋，當給定其體積與顆粒間距時，采用橢圓弧描述其表面形狀時邊界條件較多（將在3.2 節(jié)討論），故不能僅用一段橢圓弧描述。為此，這里采用兩段與大、小顆粒相對應(yīng)的橢圓弧來描述液橋的表面形狀函數(shù)y(x)：

式中：i=1、2 分別表示大、小顆粒；ai和bi分別為長、短半軸距；（0，ni）為橢圓的中心。依據(jù)式（2）和式（5）所示邊界條件，可求得式（19）中的幾何參數(shù)ni、ai和bi（i=1，2）：

為評價采用橢圓弧描述不等徑濕顆粒間液橋表面形狀的精度，可對比分析基于橢圓弧假定的不同顆粒半徑比R12的濕顆粒間液橋?qū)oung-Laplace 方程數(shù)值解的預(yù)測效果，如圖4所示，橢圓弧能較為準確地描述不同R12的濕顆粒間小體積液橋（VL?B≤1×10-3）在θ≤20°時的表面形狀。

由于接觸點橫坐標xci和最窄“頸部”半徑y(tǒng)0均隨顆粒-液橋接觸半徑y(tǒng)ci變化，故依據(jù)文獻［12］的無量綱方法引入以下兩個無量綱變量：

式中：χi為無量綱接觸點坐標；ηi為液橋最窄“頸部”半徑與顆粒-液橋接觸半徑之比。當液橋表面形狀確定時，χi和ηi分別滿足χi≥0 和0 ≤ηi≤1。

將式（2）、式（4）、式（5）和式（21）代入式（20），并在小變量yci/Ri處進行泰勒展開，則ni、ai和bi可近似表示為：

當兩個顆粒半徑相同時（R1=R2），ni、ai和bi的公式同文獻［12］；當ai =bi時，橢圓弧轉(zhuǎn)化為圓弧，即橢圓弧假定退化為圓弧假定，則χ與η的關(guān)系為：

圖4 液橋表面形狀函數(shù)y(x)的數(shù)值解及橢圓弧假定的預(yù)測結(jié)果（θ =20°且VL?B ≤1×10-3）

3.2 液橋表面形狀的補充幾何條件和封閉關(guān)系式基于橢圓弧假定描述等徑濕顆粒間液橋的表面形狀時，較圓弧假定不同之處在于需要構(gòu)建無量綱變量χ與η的封閉關(guān)系式來確定橢圓弧的自由參數(shù)［12］。當一對濕顆粒R1=R2時，可通過補充平均曲率在液橋最窄“頸部”與顆粒-液橋接觸點處取值相同的幾何條件來構(gòu)建封閉關(guān)系式。反之，由于其間液橋關(guān)于x=0 的最窄“頸部”處不對稱，則這種平均曲率相同的補充幾何條件不再適用。因此，對不等徑濕顆粒間的液橋表面形狀構(gòu)建封閉關(guān)系式時，需補充新的幾何條件。

如2.2 節(jié)所述，對式（6）所示Young-Laplace 方程積分可得到式（7）所示算子Λ[y](x)（相應(yīng)的毛細力沿軸向不變），故Young-Laplace 方程與式（7）是等效的。當控制液橋最窄“頸部”半徑y(tǒng)0和小顆粒-液橋接觸半徑y(tǒng)c2時，欲通過式（12）確定大顆粒-液橋接觸半徑y(tǒng)c1（若yc1的解存在），則需滿足Λ[y](xc1)=Λ[y](xc2)，由此可實現(xiàn)Λ[y](xc1)=Λ[y](0)=Λ[y](xc2)，進而依據(jù)Λ[y](x)在不同x處相同的原則補充幾何條件。為此，針對描述液橋表面形狀的兩段橢圓弧，依據(jù)Λ[y](x)在液橋最窄“頸部”與顆粒-液橋接觸點之間中點處相同可補充兩個幾何條件：

式中：i=2 時的幾何條件因在式（12）的推導(dǎo)過程中已考慮而得以滿足。 由式（24）可知Λ[y](xc1/2)=Λ[y](xc2/2)，則x=xc1、xc1/2 、0、xc2/2 和xc2處的毛細力均沿軸向相等，因此式（24）所示兩個補充幾何條件因Λ[y](x)為常數(shù)而考慮了Young-Laplace 方程的性質(zhì)。需注意：這兩個幾何協(xié)調(diào)條件僅考慮對液橋表面形狀函數(shù)的一階導(dǎo)數(shù)y′(x)，不同于文獻［12］中考慮對液橋表面形狀函數(shù)的二階導(dǎo)數(shù)y″(x)。

將式（22）所示參數(shù)ni、ai和bi代入式（19），同時聯(lián)立式（24）推得無量綱變量χi與ηi之間的封閉關(guān)系式：

式中：g(χi，ηi)具體表示為：

由于封閉關(guān)系式的函數(shù)對應(yīng)法則對大、小顆粒均為g，即g(χ1，η1)=0 、g(χ2，η2)=0 。為符號表述簡便，故下文統(tǒng)一采用g(χ，η)=0 的形式。

當χ≥0 、 0 ≤η≤1 且（θ≤60°） 時，g(χ=0，η)≤0 且?g(χ，η)/?χ≥0。因此，若給定η值，則式（26）所示χ的方程有唯一非負解，這個具有物理意義的唯一解可表示為χ=f(η)的函數(shù)形式。在固-液接觸角θ分別取0°和20°、顆粒半徑比R1/R2分別取4/3 和16、無量綱液橋體積分別取VL?B≤1×10-6和≤1×10-3條件下，整理χ與η關(guān)系的Young-Laplace 方程數(shù)值解與式（25）所得χ=f(η)函數(shù)及圓弧假定（式（23））的預(yù)測曲線，如圖5所示：當η=η?時，封閉關(guān)系χ=f(η)存在峰值χ?，可將(η?，χ?)稱為界限點；當χ值一定時，η=f-1(χ)存在兩段，亦如文獻［12，15］所述，當η≥η*時，表面積較小的液橋才具有物理意義。由圖5（a）（c）可知：當VL?B≤1×10-6時，大、小顆粒對應(yīng)的Young-Laplace方程數(shù)值解可壓縮為一條數(shù)據(jù)線，其與式（25）所得χ=f(η)函數(shù)的預(yù)測曲線較為吻合；大、小顆粒對應(yīng)的界限點(η*，χ*)完全相同。由圖5（b）（d）可知，當VL*B≤1×10-3時，Young-Laplace 方程數(shù)值解呈一定范圍的數(shù)據(jù)帶（并非一條數(shù)據(jù)線），大顆粒較小顆粒數(shù)據(jù)帶更窄；式（25）所得χ=f(η)函數(shù)的預(yù)測曲線較數(shù)值解偏小，Kruyt 等［12］對等徑濕顆粒間液橋亦得到了類似的計算結(jié)果。不僅如此，由圖5亦可知，當θ≤20°且VL?B≤1×10-3時，式（25）所得χ=f(η)函數(shù)較式（23）所示圓弧假定的預(yù)測曲線對不同值的Young-Laplace 方程數(shù)值解更吻合。

當計算小體積液橋（V LB≤1×10-3）的毛細力時，需要在η≥η?的穩(wěn)定段內(nèi)對η=f-1(χ)的封閉關(guān)系進行近似簡化。這里將η=f-1(χ)的封閉關(guān)系式近似為：

圖5 θ 、VL*B 和R1/R2 不同時 χ 與η 關(guān)系的數(shù)值解以及由式（25）和式（23）所得預(yù)測曲線

式中：函數(shù)L(χ)可表示為χ的二次多項式：

式中：k0、k1和k2為方程系數(shù)。

式（27）所示近似關(guān)系需滿足χ=f(η)函數(shù)的4 個性質(zhì)：①f(1)=0 ；②（由圖5可知，與式（23）所示圓弧假定在此處的斜率相同）；③f(η*)=χ*；④f′()η =0 。因此，式（27）與性質(zhì)①、②和③對應(yīng)的約束條件可表示為：

將式（30）代入式（28）所示函數(shù)L(χ)中，進而將其代入式（27）即可得到η=f-1(χ)的近似關(guān)系：

3.3 顆粒-液橋接觸半徑關(guān)于液橋體積和顆粒間距的解析公式若控制液橋體積VLB和顆粒間距2d，為推導(dǎo)適用于小體積液橋的顆粒-液橋接觸半徑（yc1和yc2）解析公式，可將式（19）所示形如橢圓弧的液橋表面形狀函數(shù)y(x)代入式（16）—（18）所示液橋體積公式，進而在yc1和yc2較小的條件下對VLB在ξ處得到泰勒展開式：

式中函數(shù)j(η)和h(η)可表示為：

在η≥η?的η=f-1(χ)穩(wěn)定段內(nèi)，j(η)和h(η)可在(χ=0，η=1) 處得到泰勒展開式：

當給定VLB和2d時，式（32）中的yc1和yc2未知，對其求解則需增加yc1與yc2的約束條件，可引入表征顆粒-液橋接觸點處毛細力相同的式（11）表示。當yc1?R1且yc2?R2時，VLB較?。╒L?B≤1×10-3），又故可在yc2和y0確定時，將式（11）簡化為關(guān)于yc1的二次方程：

當液橋體積很小時，式（11）和式（35）的解基本相同。將y0=η2yc2（由式（21）可得）代入式（35），并在小變量yc2處得到泰勒展開式：

當0≤θ≤45°且0.4≤η2≤1.0 時，函數(shù)k(η2)在0 ～2 范圍內(nèi)變化。由于yc2?R2≤R1，式（36）中的項可忽略不計，故此時大、小顆粒與液橋的接觸半徑近似相等，即yc1≈yc2，這與圖4所示小體積液橋（VL?B≤1×10-3）與大、小顆粒的接觸特征一致。由該近似條件可知，η1≈η2（依據(jù)式（21）），進而依據(jù)式（25）所示封閉關(guān)系可知χ1=f(η1)≈f(η2)=χ2。因此，可將式（32）進一步整理為VLB關(guān)于yc2的公式：

式中：Re為顆粒等效半徑，同式（1）。由yc1≈yc2和η1≈η2可知xc1≈-xc2（依據(jù)式（21））。將該近似條件代入式（4）和式（18），即可得到小顆粒至x=0 處的間距d2（圖2）公式：

將式（34）代入式（37），并聯(lián)立式（4）、式（21）、式（30）和式（38），把式（37）進一步整理為VLB和2d一定時關(guān)于yc2的五次方程：

式中：ξs為關(guān)于ξ的公式縮略代換變量：

當0≤θ≤60°且yc2?Re時， -π/2 ≤ξs≤-3/10，則式（39）中yc2的五次方項較其四次方項可忽略，而yc2的三次方項較其平方項可忽略，由此可將式（39）簡化為四次方程：

依據(jù)文獻［12］，式（41）關(guān)于yc2的解可在顆粒間距較小和較大兩個范圍內(nèi)分別求取。采用文獻［12］所述聯(lián)立方法，將式（41）的解表示為：

式中：μ（不同于文獻［12］）、?、w1(μ)和w2(?)分別表示為：

需注意：如式（42）所示，yc1和yc2僅通過Re反映了大、小顆粒半徑（R1和R2）的影響。

3.4 小體積液橋斷裂距離的解析公式已有試驗［10，16］表明：濕顆粒間的液橋會隨顆粒間距增大而逐漸拉伸直至斷裂，而液橋斷裂距離2dr與其體積VLB和固-液接觸角θ有關(guān)，并與Young-Laplace 方程出現(xiàn)唯一數(shù)值解時的數(shù)據(jù)點對應(yīng)［9-10，12］。依據(jù)文獻［12］， 2dr隨其VLB和θ的變化關(guān)系亦可采用橢圓弧假定預(yù)測。

由圖5中式（25）所示封閉關(guān)系χ=f(η)可知，χ在0≤η≤1 范圍內(nèi)存在峰值，該峰值點（或界限點）的坐標用(η*，χ*)表示，相應(yīng)的j(η)值和h(η)值（函數(shù)公式如式（33）所示）分別用j?和h?表示。由于不等徑濕顆粒間液橋即將斷裂時其體積很小，yc2?Re，則式（37）中yc2的四次方項可忽略，由此可將式（37）所示VLB的公式進一步簡化為：

此外，當(η，χ)位于界限點(η?，χ?)時，yc較小且d=dr，聯(lián)立式（4）（21）（38）可得χ?關(guān)于dr和yc2的近似公式：

由此，聯(lián)立式（44）和式（45）可得到2dr的解析公式：

式中：不同θ值對應(yīng)的χ?和η??值（由式（25）確定）、j?和h?值（由式（33）確定）如表1所示。此外，對比文獻［9］和式（46）中系數(shù)α在0≤θ≤20°范圍內(nèi)的取值，如圖6所示，式（46）計算的α值與文獻［9］差異較小，不超過7%。

3.5 小體積液橋毛細力的解析公式針對不等徑濕顆粒間形成的小體積液橋，采用橢圓弧假定，推導(dǎo)液橋毛細力關(guān)于其體積和顆粒間距的解析公式（無需引入任何標定的擬合參數(shù)）。由式（10）可知，計算毛細力F的關(guān)鍵在于求解毛細長度系數(shù)λ，λ可通過液橋最窄“頸部”半徑y(tǒng)0、顆粒-液橋接觸半徑（如yc2）及其接觸點斜率（如y′c2）計算得到。若控制液橋體積VLB、顆粒間距2d和固-液接觸角θ，當VLB較小時，顆粒-液橋接觸半徑（yc1≈yc2）可由式（42）確定；小顆粒至x=0 處的間距d2（圖2）由式（38）確定，再通過式（4）和式（21）確定χ2，進而將χ2、χ?和η?值（依據(jù)表1查取的）代入式（27）以確定η2；將η2代入式（21）即可確定y0，繼而由式（5）確定y′c2；將y0、yc2和y′c2代入式（8），并在yc2處進行泰勒展開，從而得到λ的解析公式：

最后，將式（47）代入式（10）確定F；若采用顆粒等效半徑Re對F無量綱化，結(jié)合式（47）亦可確定無量綱毛細力F*（F=λ/Re）。

整理不同顆粒半徑比（R1/R2=2 和128）的濕顆粒間液橋取不同無量綱體積（體積較小，VL?B=1×10-6和1×10-3）和顆粒間距之比d?/dr?（d*和dr?分別為無量綱的顆粒半間距和液橋斷裂半倍距離）時無量綱毛細力F*的Young-Laplace 方程數(shù)值解及式（47）的預(yù)測值，如圖7所示（R1/R2=128 可近似視為球體與片狀濕顆粒之間形成的液橋），式（47）對Young-Laplace 方程數(shù)值解的預(yù)測較為吻合：當d*＜0.6dr?時，VL?B=1×10-6和1×10-3時的預(yù)測誤差分別小于5%和6%；當d*≥0.6dr?時，預(yù)測誤差雖有所增大，但此時F*已很小。原因可能在于：當顆粒間距較大時，大、小顆粒與液橋的接觸半徑差異顯著，因此依據(jù)yc1≈yc2條件所計算的F*值較Young-Laplace 方程數(shù)值解勢必會產(chǎn)生誤差。此外，當d*≥0.6dr?時，極小體積液橋（VL?B=1×10-6）對應(yīng)的Young-Laplace 方程數(shù)值解及式（47）的預(yù)測值會出現(xiàn)波動（4.2 節(jié)的圖12（a）亦有此現(xiàn)象）。原因可能在于：液橋表面的自由能對其穩(wěn)定性影響顯著［9］，當極小體積液橋被拉伸接近斷裂時，會因其表面（尤其在最窄“頸部”處）的自由能過小而無法控制其穩(wěn)定性，致使波動現(xiàn)象出現(xiàn)；當液橋體積增大時，其表面自由能亦隨之增大，從而使這種波動現(xiàn)象逐漸消失（4.2 節(jié)的圖12（b）—（d）可佐證）。

圖6 文獻［9］和式（46）中系數(shù)α 與θ 的關(guān)系

表1 不同θ 值對應(yīng)的 χ *、η*、 j *和h*值

不僅如此，為從試驗角度研究濕顆粒間不同體積液橋（液橋體積亦較小，VL?B≤1×10-3）的毛細力隨顆粒間距的變化特征，Rabinovich 等［16］利用原子力顯微鏡（AFM）測定了3 對玻璃質(zhì)球體濕顆粒被液體不完全浸潤（θ=10°）時其間形成的液橋（液體選取石蠟油，其表面張力為27.5×10-3N/m）在不同顆粒間距對應(yīng)的毛細力，其中R1/R2≈1.84、1.71和1.44；VL?B≈1.3×10-5、8.7×10-5和3.2×10-4；相應(yīng)的試驗編號分別為①、②和③。整理這3 對不等徑球體濕顆粒間液橋取不同d*/dr*時F*的實測值及式（47）的預(yù)測值如圖8所示，式（47）對試驗①、②、③中液橋F*的預(yù)測值與其實測值較為吻合，預(yù)測誤差不超過9%，僅當d*≥0.5dr?時，預(yù)測誤差雖有所增大，但此時F*已很小，與圖7所示預(yù)測結(jié)果類似。

圖7 R1/R2、和d */dr* 不同時無量綱毛細力F*的Young-Laplac e 方程數(shù)值解及式（47）的預(yù)測值（θ =20°）

圖8 和不同時F *的實測值及式（47）的預(yù)測值（θ =10°）

當給定VLB、2d、θ和顆粒半徑（R1和R2）時，則可依據(jù)橢圓弧假定按照如下方法計算不等徑濕顆粒間小體積液橋（VL?B≤1×10-3）的F：①檢查2d是否小于式（46）所確定的2dr；②依據(jù)式（42）由確定的VLB和d計算yc2；③將yc2代入式（38）和式（4），分別計算得d2和xc2；④采用式（21）和式（27）分別得到χ2和η2；⑤聯(lián)立式（47）和式（10）即可計算出F。

如3.3 節(jié)依據(jù)橢圓弧假定的分析所述，yc1≈yc2，進而據(jù)式（37）所示液橋體積公式可知，yc1和yc2通過式（1）所示Re僅與R1、R2有關(guān)。由此可見，當液橋體積較小時，針對等徑濕顆粒推得的液橋斷裂距離與其毛細力的解析公式可通過Re推廣至不等徑濕顆粒的情況。這種思路雖已通過對Young-Laplace 方程數(shù)值解的分析得以應(yīng)用［10-11，30］，但鮮從解析角度得到論證，為此，本節(jié)在橢圓弧假定的小體積液橋適用范圍內(nèi)從解析角度驗證了這種思路的準確性，而在大體積液橋斷裂距離與其毛細力公式中采用Re的適用性將在下節(jié)闡述。

4 基于擬合公式計算大體積液橋的斷裂距離及其毛細力

不等徑濕顆粒間大體積液橋的斷裂距離與其毛細力研究報道較少［10-11］，Willett 等［10］選取有限的Young-Laplace 方程數(shù)值解數(shù)據(jù)組構(gòu)建了液橋斷裂距離及其毛細力的擬合公式，Lian 等［11］基于任意顆粒間距時的毛細力與顆粒相互接觸（顆粒間距為零）時的初始毛細力之比構(gòu)建了液橋毛細力擬合公式。然而，由于相互接觸的濕顆粒間無法形成較大體積的凹形鐘擺狀液橋，因此他們［10-11］提出的公式僅適用于小體積液橋，對大體積液橋適用性有限。為此，本節(jié)以大體積液橋為研究對象，旨在分析顆粒半徑比和固-液接觸角對其斷裂距離與毛細力的影響。為對比文獻［10-11］的研究，本節(jié)采用較文獻［10］規(guī)模更大的Young-Laplace 方程數(shù)值解數(shù)據(jù)組，選取較文獻［11］范圍更廣的液橋體積和顆粒半徑比。

當液橋體積和顆粒間距較大時，橢圓弧假定不適于描述其表面形狀。為此，基于2.3 節(jié)所述Young-Laplace 方程數(shù)值解的數(shù)據(jù)組，采用顆粒半徑比而非顆粒等效半徑來描述顆粒半徑的影響，一方面構(gòu)建液橋斷裂距離的擬合公式，另一方面將Lian 等［11］提出的毛細力公式進行改進，使之適用于大體積液橋。

4.1 液橋斷裂距離的擬合公式：從小體積液橋到大體積液橋Willett 等［10］從試驗與數(shù)值解角度分析了顆粒半徑比對液橋斷裂距離的影響，依據(jù)試驗結(jié)果闡述了顆粒半徑比對大體積液橋斷裂距離的影響，進而對有限的Young-Laplace 方程數(shù)值解數(shù)據(jù)組進行曲線擬合，得到了無量綱液橋斷裂距離2dr?的擬合公式：

式中：固-液接觸角θ用弧度表示；VL?B為無量綱液橋體積，其與2dr?的無量綱定義方式如2.3 節(jié)所述。對比分析不同顆粒半徑比（1≤R1/R2≤128）的濕顆粒間不同體積（0≤≤0.22）液橋在θ=20°、40°時2dr?的Young-Laplace 方程數(shù)值解及式（48）的預(yù)測值，如圖9（a）（c）所示，式（48）對θ=20°、40°時數(shù)值解的預(yù)測誤差較大（分別為4.5%和10%）。因此，有必要在式（48）基礎(chǔ)上通過引入和θ的高階項以構(gòu)建適用液橋體積更廣范圍的無量綱液橋斷裂距離擬合公式。通過對2dr?的Young-Laplace 方程數(shù)值解進行曲線擬合可得到2dr?關(guān)于VL?B、θ和R1/R2的擬合公式：

對比分析式（49）對不同顆粒半徑比（1≤R1/R2≤128）的濕顆粒間不同體積（0≤VL*B≤0.22）液橋在θ=20°、40°時2dr?的Young-Laplace 方程數(shù)值解的預(yù)測結(jié)果，如圖9（b）（d）所示，式（49）對θ=20°、40°時數(shù)值解的預(yù)測誤差較式（48）顯著減小，分別為1%和2%。

圖9 采用式（48）與式（49）對不同R1/R2 和θ 時2dr? 的Young-Laplace 方程數(shù)值解的預(yù)測結(jié)果比較

4.2 液橋毛細力的擬合公式評價：從小體積液橋到大體積液橋Willett 等［10］和Lian 等［11］已提出大體積液橋的毛細力擬合公式，這里采用本文構(gòu)建的Young-Laplace 方程數(shù)值解數(shù)據(jù)組評價其預(yù)測精度。

Willett 等［10］通過對Young-Laplace 方程數(shù)值解進行曲線擬合，構(gòu)建了無量綱液橋毛細力F*的擬合公式（F*=λ/Re）：

Lian 等［11］先分析了濕顆粒相互接觸（顆粒間距為零）時的初始毛細力，進而依據(jù)任意顆粒間距與零間距時的毛細力比值構(gòu)建了F*的擬合公式：

式中dθ可表示為：

式中：2dr?可通過式（48）令R1/R2=1 而確定；β0為充填角（用弧度表示），可采用給定VL?B和θ時關(guān)于β0的隱式公式確定：

依據(jù)大體積液橋F*的Young-Laplace 方程數(shù)值解（如圖10所示，VL?B=0.13 且θ=20°）可知：F*（已采用顆粒等效半徑無量綱化）與VL?B、 2d?、和θ有關(guān)，即如果顆粒半徑（R1和R2）的影響可用Re描述，則F*應(yīng)與無關(guān)。當2d?較小時，對F*基本無影響（圖10），這與文獻［10］中的研究結(jié)論一致；當2d?較大甚至液橋即將斷裂時，對F*影響較為顯著（圖10）。進一步整理不同條件下液橋恰好斷裂時Fr?（即F?(2d=2dr?)）與VL?B的關(guān)系，如圖11所示：當液橋體積較小時（即VL?B≤1×10-3），對Fr?基本無影響；當液橋體積較大時（即VL?B＞1×10-3），對Fr?影響顯著。因此，當液橋體積較小時，毛細力公式中無需考慮的影響；反之，當液橋體積較大時，則應(yīng)在該公式中考慮的影響。

圖10 R1/R2 不同時F *與2d *的關(guān)系

圖11 R1/R2 不同時Fr *與VL*B 的關(guān)系

由此可見，式（52）（53）所示無量綱的液橋毛細力及其斷裂距離擬合公式均未考慮的影響，因而Lian 等［11］的擬合公式無法精確描述不等徑濕顆粒間大體積液橋的毛細力。為對比分析Wil?lett 等［10］和Lian 等［11］提出的毛細力擬合公式對Young-Laplace 方程數(shù)值解的預(yù)測精度，選取R1/R2=16 且θ=20°時，不同體積（VL?B=1×10-6、1×10-3、0.03 和0.09）液橋的F*與關(guān)系的Young-La?place 方程數(shù)值解及這兩種公式預(yù)測值作為示例，如圖12所示。由圖12（b）—（d）可知，當d*≤0.5dr?、θ≤40°且1×10-3≤VL?B≤0.13 時，Willett 公式的預(yù)測誤差小于5%；然而，對比圖12（a）（b）可知，當VL?B≤1×10-3時，其預(yù)測誤差超過10%。相比之下，當d*≤0.5dr?且VL?B≤1×10-3時，Lian 公式的預(yù)測誤差小于10%。然而，這種公式對大體積液橋預(yù)測誤差較大，由圖12（c）（d）可知，當VL?B=0.09 時，Lian 公式因其未考慮R1/R2的影響而對F*的Young-Laplace 方程數(shù)值解預(yù)測誤差超過10%。此外，由圖12（a）亦可看出：極小體積液橋（VL?B=1×10-6）在臨近斷裂的過程中，相應(yīng)的Young-Laplace方程數(shù)值解及不同種公式的預(yù)測值會出現(xiàn)不同程度的波動現(xiàn)象，這與3.5 節(jié)所述液橋表面（尤其在最窄“頸部”處）的自由能對其穩(wěn)定性的顯著影響有關(guān)，而且當液橋表面的自由能隨其體積增大時，這種波動現(xiàn)象會消失（如圖12（b）—（d）所示）。

圖12 VL*B 和不同時F *的Young-Laplace方程數(shù)值解與Willett公式、Lian公式及其引入式（49）的改進公式預(yù)測值（=16和θ =20°）

為提高Lian 公式的預(yù)測精度，可在式（52）所示F*公式或式（53）所示dθ公式中引入R1/R2。出于計算合理、簡便考慮，采用預(yù)測無量綱液橋斷裂距離更為合理的式（49）替換式（53）所示2dr?，以考慮R1/R2的影響。由圖12可知，改進的Li an 公式對小體積液橋至大體積液橋（≤0.13）的毛細力Young-Laplace 方程數(shù)值解預(yù)測均較為吻合，預(yù)測誤差小于5%，優(yōu)于兩種已有公式。由圖12（a）亦可知，不等徑濕顆粒間的大體積凹形鐘擺狀液橋僅存于最小顆粒間距大于零范圍內(nèi)。當給定小顆粒-液橋接觸半徑y(tǒng)c2和液橋最窄“頸部”半徑y(tǒng)0時，該最小顆粒間距對應(yīng)2.2 節(jié)所述大顆粒-液橋接觸半徑y(tǒng)c1解存在區(qū)域的下邊界。由此可見，文獻［11］所述采用零顆粒間距進行變量無量綱化的方法不適用于大體積液橋（如VL?B= 0.09）。

4.3 考慮顆粒粒徑和液橋體積的毛細力計算方法應(yīng)用討論本構(gòu)模型的構(gòu)建是土壩變形及其穩(wěn)定性計算中的關(guān)鍵環(huán)節(jié)。如研究背景所述，當壩前水位驟降時，構(gòu)成壩體的濕顆粒材料內(nèi)部的毛細黏聚作用至關(guān)重要。如文獻［8］所述，要在濕顆粒材料的微細觀結(jié)構(gòu)彈塑性本構(gòu)模型中考慮這種毛細黏聚作用，則需計算濕顆粒間的毛細力，文獻［8］中雖引入了毛細力公式，但僅與顆粒間距有關(guān)，因未考慮不同體積液橋而無法考慮水位驟降引起飽和度變化對此類本構(gòu)模型的影響。為此，若以顆粒級配較分散的濕顆粒材料為研究對象，將其簡化為不等徑濕顆粒-液橋模型，顆粒半徑（R1和R2）分別取這種材料最大、最小粒徑的一半，固液接觸角θ≤40°，再采用等效顆粒半徑Re對液橋體積VLB和顆粒間距2d進行無量綱化得到相應(yīng)的無量綱變量VL?B和2d?，同時把液橋劃分為小體積（VL?B≤1×10-3）和大體積（1×10-3＜VL?B≤0.13）兩個范圍，則可將本文提出的考慮顆粒粒徑和大、小體積液橋的毛細力計算公式嵌入微細觀結(jié)構(gòu)彈塑性本構(gòu)模型中，具體可依據(jù)圖13所述思路從以下4 個步驟實現(xiàn)：

圖13 不同體積液橋斷裂距離及其毛細力計算公式在微細觀結(jié)構(gòu)彈塑性本構(gòu)模型中的應(yīng)用思路

（1）當VL?B≤1×10-（3θ≤20°）時，依據(jù)3.5 節(jié)所述將液橋斷裂距離解析公式（式（46））與毛細力解析公式（式（47））相結(jié)合的思路，計算小體積液橋產(chǎn)生的粒間毛細力F。

（2）當1×10-3＜VL?B≤0.13 時，依據(jù)4.2 節(jié)所述引入式（49）改進Lian 等［11］的毛細力擬合公式（式（52））的思路，計算大體積液橋產(chǎn)生的粒間毛細力F。

（3）計算不同體積液橋在濕顆粒構(gòu)成孔隙空間內(nèi)所占的飽和度：依據(jù)圖2所示不等徑濕顆粒之間構(gòu)成的孔隙空間計算孔隙體積Vv，則液橋所占的飽和度可按照下式計算：

（4）將粒間毛細力嵌入微細觀結(jié)構(gòu)彈塑性本構(gòu)模型進行土壩變形及其穩(wěn)定性計算的基本思路：將第（1）、（2）步所述不同體積液橋所產(chǎn)生的粒間毛細力作為外荷載加入濕顆粒間的接觸力中，進而將所有粒間接觸力代入文獻［8］所述Hertz-Mindlin 彈性接觸公式和Mohr-Coulomb 屈服準則中，以分析第（3）步所計算飽和度對這種材料毛細黏聚作用的影響，進而結(jié)合非飽和土的水力特性參數(shù)（持水曲線、滲透系數(shù)函數(shù)）與有效應(yīng)力張量在離散單元（DEM）程序中建立土壩在非飽和滲流作用下的數(shù)值模型，即可開展土壩在水位驟降時的變形及其穩(wěn)定性計算。這也是本文下一步的研究方向。

5 結(jié)論

為研究濕顆粒間的細觀毛細黏聚作用機理，開展了考慮顆粒粒徑和液橋體積的毛細力計算分析，將濕顆粒簡化為一對不等徑球體顆粒，結(jié)合其間液橋的幾何特征構(gòu)建了Young-Laplace 方程數(shù)值解的數(shù)據(jù)組，一方面基于橢圓弧假定對小體積液橋推導(dǎo)解析公式，另一方面基于Young-Laplace 方程數(shù)值解構(gòu)建適用于大體積液橋的擬合公式，主要得出以下結(jié)論：

（1）橢圓弧假定對不等徑濕顆粒間的小體積液橋（VL?B≤1×10-3）在固-液接觸角θ≤20°的范圍內(nèi)是適用的。由該假定可知，大、小顆粒與小體積液橋的接觸半徑基本相同（即yc1≈yc2），并且已通過Young-Laplace 方程數(shù)值解得到驗證。當液橋體積較小時，依據(jù)橢圓弧假定的解析研究和Young-Laplace 方程數(shù)值解分析表明：采用顆粒等效半徑可將等徑濕顆粒的解析公式推廣至不等徑濕顆粒的范圍?；跈E圓弧假定推得液橋斷裂距離及其毛細力關(guān)于液橋體積和顆粒間距的解析公式，并采用已有文獻中不等徑球體顆粒間小體積液橋的毛細力實測值驗證了解析公式的有效性。

（2）當液橋體積較大時（VL?B＞1×10-3），由Young-Laplace 方程數(shù)值解可知：不等徑濕顆粒的公式不能簡單地通過在等徑濕顆粒公式中采用顆粒等效半徑代換得到。因此，對液橋體積變化更廣范圍內(nèi)Young-Laplace 方程數(shù)值解的數(shù)據(jù)組進行曲線擬合，通過引入顆粒半徑比構(gòu)建了適用于大體積液橋（VL?B＞1×10-3）的斷裂距離擬合公式。

（3）將（2）提出的斷裂距離擬合公式引入已有不等徑濕顆粒間液橋的毛細力擬合公式進行改進，并采用Young-Laplace 方程數(shù)值解的數(shù)據(jù)組評價了已有毛細力擬合公式的適用性，通過對比預(yù)測誤差發(fā)現(xiàn)：改進公式對半徑比在1～128 范圍內(nèi)的濕顆粒間不同體積液橋（VL?B≤0.13）在固-液接觸角θ≤40°且顆粒間距不超過液橋斷裂距離范圍內(nèi)的毛細力預(yù)測效果優(yōu)于已有公式。

本文提出的毛細力計算方法可嵌入典型濕顆粒材料的微細觀彈塑性本構(gòu)模型，旨在用于水位驟降時土壩在非飽和滲流作用下的變形分析及其穩(wěn)定性評價。