陸培毅,韓亞飛,王成華

(天津大學 建筑工程學院，天津 300350)

0 引 言

土體流固耦合是研究土體滲流場和應力場之間的相互作用，涉及到土力學，流體力學和兩學科的交叉。太沙基最早提出有效應力原理并基于該原理提出了一維固結理論。此后，比奧基于彈性理論，平衡條件，連續(xù)條件等提出了真三維固結理論[1]，以上兩種理論是土力學中最具有代表性的流固耦合理論。后來，諸多學者分別從耦合形式，有限元等方面對以上兩種理論進行了改進和應用，形成了現(xiàn)在的流固耦合理論。假設條件和耦合角度的不同，得出耦合理論的適用范圍和精度也不盡相同。因此有必要弄清各種耦合理論的前提條件、內在聯(lián)系、區(qū)別和優(yōu)缺點。所以筆者對現(xiàn)有的耦合理論從計算域耦合、基本未知數(shù)耦合、參數(shù)耦合和本構關系耦合4個方面進行分類，目的在深入了解每種理論并總結優(yōu)缺點，以便在理論指導實際工程時做到揚長避短。

1 計算域耦合

1.1 概 述

計算域之間的耦合是指固體域和流體域之間的相互作用，按作用機理可以分為兩種，一種是固體域和流體域部分或完全重疊在一起。該種問題的求解需要分別建立固體變形和流體壓力的方程，然后根據(jù)兩者之間的變形協(xié)調或者兩者之間的相互作用力建立耦合關系。Terzaghi一維固結理論就是基于這種方法建立的。另一種是固體域和流體域完全分開，既耦合作用僅發(fā)生在兩相交界面處。這種問題的求解需要建立起兩相耦合面上力的平衡和位移協(xié)調的關系。問延煦等[2]提出的附加質量法以及水體位移元法就是基于這種方法建立的。

1.2 理論發(fā)展

太沙基一維固結理論是在一系列假設上建立起來的，倫杜利克把太沙基的一維固結理論推廣到二維和三維，后來的研究者對太沙基一維固結理論進行改進、推廣，使其更符合工程實際。改進一般基于土體的各向異性、非達西滲流和土的彈塑性本構等。問延煦等[2]已發(fā)表了太沙基一維固結理論的研究綜述，且內容詳實，此處不再贅述。

一維固結理論的不足已經被人們熟知，例如假設土體彈性，滲透系數(shù)不變等，但由于其原理清晰，計算簡便，一直被應用至今。要深刻理解其應用假設和前提，以免造成工程事故。

羅曉輝[3]對滲流場進行了穩(wěn)定滲流與非穩(wěn)定滲流有限元分析，將滲流力的作用等效為節(jié)點荷載施加到劃分好的土體上進行分析，然后將應力狀態(tài)的變化帶入滲透系數(shù)與應力狀態(tài)的關系式，從而實現(xiàn)應力場和滲流場的耦合。

使用有限元方法，將基坑降水引起的滲流力等效成節(jié)點荷載施加到土體單元網格上來模擬滲透力對土體有效應力的影響，反之，將土中有效應力的變化與滲透系數(shù)建立起聯(lián)系，實現(xiàn)了流固耦合分析，并以此來分析大面積降水引起的基坑穩(wěn)定問題。

將滲流力簡化為體力施加到土體有限元節(jié)點方法的物理意義較為明確且符合實際，將其與土體應變-滲透系數(shù)之間的關系聯(lián)合應用就可以完整的模擬流固耦合現(xiàn)象，有較高的應用價值。較為困難的是根據(jù)不同的土質建立起有效應力和滲透系數(shù)之間的關系，有待進一步研究。

H. M. WESTERGAARD[4]提出了動水壓力附加質量法公式，他將壩面視為直立且剛性的，大大化簡了水與壩體之間的流固耦合作用。該理論存在假定，只是壩體和水的一種理想化模型，但反應了動水壓力的一些本質特征且有計算簡便的優(yōu)點，適合對精度要求不高的小型工程進行初步估算。

(1)

式中：Pw(h)為水深h處的動水壓力;ah為水平設計地震加速度代表值;ρ為水體質量密度;H0為庫水深。

CLOUGH推廣了附加質量法公式，使其適用于任意形狀的壩體，并提出了適用于有限元分析的廣義質量表達式：

(2)

式中：Mai為附加質量；λi是壩面點的法向矢量；Ai為該點在壩面上的隸屬面積。

P. CHAKRABARTI等[5-6]分別在1972年和1973年對附加質量模型做出了修正，改進后的模型可以反應壩面的彈性，這樣就更接近壩水流固耦合實際。

朱驍健等[7]在動水附加質量法公式的基礎上，考慮庫水的輻射阻尼耗能，采用庫水阻尼元件法將庫水慣性質量和阻尼元件共同作用在大壩迎水面上。阻尼元件的阻尼系數(shù)根據(jù)能量在邊界的反射與透射得到，由庫水密度、水中縱波速度、單個阻尼元件的影響面積和邊界輸入、輸出能量的比值決定。將結果與實際對比，考慮水的可壓縮性更符合工程實際。

動水壓力法是建立在諸多假設上得出的，例如采用剛性壩面假設，忽略了壩體的彈性變形和庫水的狀態(tài)對動力響應的影響，這就導致其不是真正的流固耦合，且一般認為動水附加質量法夸大了動水壓力作用，與實際情況存在較大差距。

E.L.WILSON等[8]首次使用水體位移元法來求解水和壩體之間的流固耦合問題，與動水壓力法相比，該種方法將庫水和壩體都以結點位移為主變量，使其統(tǒng)一在一個系統(tǒng)中求解壩體與庫水的動力耦合問題。

陳和群等[9]通過限制水體的旋轉變量實現(xiàn)無旋性，用水體位移元法，建立壩體與庫水動力耦合的有限元法分析模型。通過分析二灘雙曲拱壩與庫水的耦合作用后發(fā)現(xiàn)水體位移元法能較好地模擬壩與庫水動力實際的相互作用。

水體位移元法同時將水的節(jié)點位移作為主變量，因此考慮了水體的可壓縮性，但目前研究較少,可以進行深入研究。

在計算域之間耦合的各種計算方法大致可以分為滿足變形協(xié)調條件的太沙基一維固結理論、將滲流力簡化為體力施加到土體有限元節(jié)點的方法、動水壓力法和水體位移元法。其中，滲流力簡化為節(jié)點力的方法可以結合下文中的基本未知數(shù)耦合理論進一步完善。水體位移法也應進一步研究。

2 基本未知數(shù)耦合

2.1 概 述

基本未知數(shù)之間的耦合指事先假定若干個未知量，然后根據(jù)耦合現(xiàn)象必須滿足的物理力學條件，邊界條件等列出方程，求解出未知量，從而達到耦合的目的。此種方法以比奧固結理論為代表，后人在其理論基礎上進行完善。SAVAGE和BRADDOCK將三維比奧固結理論應用于對各向同性的孔隙彈性介質的分析。ZIENKIEWIEZ和SHIOMI在比奧固結理論的基礎上，考慮了多孔介質的幾何、邊界及材料的非線性?？傊?，對比奧固結理論的完善主要是從有限元、本構模型、非達西滲流三方面開展的。

2.2 理論發(fā)展

比奧固結理論[1]先假定4個未知量(3個方向的位移和超靜孔隙水壓力)，然后通過彈性理論，線彈性應力應變關系，水流連續(xù)方程解出未知量，從而得到應力場和滲流場之間的耦合關系。

R. S. SANDHU等[10]最早使用有限元法來求解比奧固結方程。

在國內，沈珠江應用變分原理首先把比奧固結理論的有限單元法應用于固結分析。隨后，鄧岳保等[11]分別采用不同的方法推導了比奧固結理論有限元方程。

比奧固結理論假定土體為線彈性體，平揚等[12]將土體本構關系中的彈性矩陣替換為彈塑性矩陣，這就將該理論推廣到彈塑性分析領域，并通過水土特征曲線求出自由面的滲透系數(shù)，改進了比奧固結理論中不考慮非飽和效應的假設。

王媛[13]以比奧理論為基拙 , 假定土體飽和且土骨架滿足線彈性各向同性體。使用伽遼金法對比奧固結方程進行了求解,同時考慮了正交各向異性的達西滲流。

對于比奧固結理論的改進，無論是將土體視為各向同性還是各向異性，其分析原理都是研究外荷載下孔隙水壓力與有效應力(或總應力)及相應的變形之間的關系。因而，這種分析方法其實與滲透力作用無關，所探討的滲流場實質是孔隙水壓力分布場，并不是滲流作用下所形成的滲透力場。

固結理論的非達西滲流多集中在一維固結理論的研究中，而對于比奧固結的非達西滲流研究較少。太沙基在1925年指出，達西定律不適用于塑性大的黏性土。劉慈群[14]在PASCAL的基礎上，得出了上述問題的近似解。劉忠玉等[15]引入可以同時考慮低速滲流曲線段和較高速滲流直線段的非達西滲流方程，重新推導飽和黏土一維固結方程，并采用有限體積法對該方程進行數(shù)值求解。并探討非達西滲流參數(shù)對固結過程的影響。計算結果表明，非達西滲流延緩了飽和黏土中孔隙水壓力的消散速度，故使得地基的固結速度比太沙基一維固結理論值要慢，最后討論了Terzaghi一維固結理論的適用范圍。李傳勛等[16]研究了基于指數(shù)形式、非牛頓指數(shù)的一維固結理論。

S. HANSBO[17]最早對非達西滲流進行了研究，他首次推導了考慮非達西滲流的豎井地基固結解析解。S. HANSBO的非達西滲流模型表達式：

(3)

式中：i0為起始水力梯度；iL為門檻水力梯度。

T. C. ING等[18]基于虛功原理推導了HANSBO非達西滲流模型的軸對稱比奧固結有限元方程，并分析了非達西滲流對固結計算的影響。除了得出非達西滲流對固結的延緩作用這一普遍結論，還得出，只有在iL≥40且m≥1.5時，非達西滲流才對固結有明顯的影響。另外，應力歷史對固結速率重要影響。

在比奧固結理論基礎上的非達西滲流研究中，非達西滲流能更好的描述黏性土和小荷載情況下的滲流，但由于土體的成層性導致的水平向和豎向滲透系數(shù)的差異是非達西滲流不能描述的。

3 參數(shù)耦合

3.1 概述

流固耦合是一個動態(tài)過程，在實際的滲流過程中, 由于孔隙流體壓力的變化, 一方面要引起多孔介質骨架有效應力變化, 由此導致多孔介質滲透率和孔隙率的變化;另一方面, 這些變化又反過來影響孔隙流體的流動和壓力的分布。參數(shù)之間的耦合就是建立起土體應力或應變與滲透系數(shù)或孔隙比之間的關系，達到耦合的目的。其中滲透系數(shù)則是應力場、滲流場相互耦合的“橋梁”，也是實現(xiàn)真正的滲流耦合分析的關鍵。對于假定滲透系數(shù)為常數(shù)的滲流耦合分析，只能體現(xiàn)滲流場對應力場的影響，不能體現(xiàn)應力場對滲流場的影響。

3.2 理論發(fā)展

Kozeny-Craman公式建立了土體滲透系數(shù)與孔隙率的關系式：

(4)

式中：K0，n0分別為初始滲透系數(shù)和孔隙率；K，n分別為耦合分析滲透系數(shù)和孔隙率。

駱祖江等[19]在耦合計算分析時，將n=(n0+εv)/(1+εv)帶入式(4)得到滲透系數(shù)和體應變的關系。

李培超等[21]基于考慮孔隙率和孔隙流體壓力的多孔介質有效應力原理，推導出孔隙率和滲透系數(shù)的動態(tài)變化關系，然后建立應力場和滲流場方程，得到比較完善的流固耦合數(shù)學模型。

在滲透系數(shù)和孔隙比的關系中，滲透系數(shù)是隨著孔隙率變化的，因此用統(tǒng)一的滲透系數(shù)來表示整個滲流場的滲流特性是不合理的。

滲透系數(shù)與孔隙比關系式匯總如表1：

表1 滲透系數(shù)與孔隙比關系匯總Table 1 Summary of relationships between void ratios and permeability coefficients

注：其中Deff為土平均有效粒徑；ek與當前應力狀態(tài)相關的空隙比；CF為土體的粘粒含量；Ac為土體活性指數(shù)；Ck為土體滲透系數(shù)相關的指標，一般Ck=0.5e0。

冉啟全等[26]研究了油藏數(shù)值模擬中的流固耦合，其中通過Kozeny-Craman方程，推導出等溫滲流過程中滲透率與體積應變的關系，建立起油藏開采過程中，油藏壓力變化和巖體變形之間的固耦合關系。

(5)

式中：K0，K意義同上；φ0為孔隙度；εv為體積應變。

陳曉平等[27]給出了非均質土壩滲流場和應力場耦合的數(shù)學模型，其方法是將滲流場和應力場分開求解，分別寫出與應力場有關的滲流場方程和與滲流場有關的應力場方程，并采用滲透系數(shù)和應變關系的經驗公式：

(6)

李筱艷等[28]根據(jù)抽水試驗及沉降觀測資料，建立了土體滲透系數(shù)與有效應力增量的非線性耦合模型，通過線性回歸得到了它們的關系式：

K=K0exp(-λ×Δσ)

(7)

式中：λ為試驗參數(shù)，反應了土體中滲透系數(shù)隨有效應力變化的幅度。

馬少坤等[29]采用ABAQUS中修正劍橋模型，推導出K0固結狀態(tài)下土體的孔隙比隨深度變化的關系，并建立滲透系數(shù)隨深度線性變化的關系，然后導入修正劍橋模型進行降水開挖流固耦合分析，并與不考慮降水的開挖分析作對比，從而得到考慮孔隙比和滲透系數(shù)隨深度變化時基坑降水開挖流固耦合作用下的圍護結構的變形及土體變形特征。

在滲透系數(shù)和體應變的關系中，對不同的應力路徑，用統(tǒng)一的關系式表示滲流系數(shù)與體應變的關系也是不符合實際的。因為在不同應力路徑條件下，即使各主應變不同，體應變也可能會相同，但是不同的主應變對滲透系數(shù)的影響不同。例如在豎向加載時，主要表現(xiàn)為ε1的變化，此時對豎向滲透系數(shù)影響較大；而在側向卸載時，主要表現(xiàn)為ε3的變化，此時對水平向滲透系數(shù)影響較大，但兩種情況下的體應變可能是一樣的。

4 本構關系耦合

4.1 概 述

土體液化和管涌現(xiàn)象均涉及流固耦合。土體一般在動力荷載作用下發(fā)生液化，研究液化機理必須研究土體中孔隙水壓力的變化以及土體和流體之間的相互作用；管涌現(xiàn)象多發(fā)生在顆粒級配不良的粉土中，土體在滲流作用下，逐漸被侵蝕轉化為液化細顆粒，液化細顆粒隨水流失導致孔隙率增大，從而增強了土體的滲透性；滲透性的增強使得滲流速度增大，滲流速度的增大進一步加劇了土體的侵蝕。即侵蝕與滲流之間存在著耦合效應，兩者相互促進，相互影響，聯(lián)系兩者的紐帶就是孔隙率。在研究該問題時需要提出具體的本構模型來反應耦合關系。

4.2 管涌多相耦合理論的發(fā)展

大量試驗和現(xiàn)場資料表明，管涌的發(fā)生是一個由骨架相、液相和液化細顆粒相之間的相互作用的過程。管涌產生機理如圖1，管涌多相耦合示意圖如圖2[30]。多相對管涌的研究經歷了不考慮流固耦合到考慮流固耦合再到考慮侵蝕本構關系耦合3個階段，使對管涌的產生機理有了更深刻的認識。

周恩全[31]研究了飽和砂土液化后的流體本構模型，及將液化后的土體當作流體來研究，是一種較新穎的求解思路，他通過實驗發(fā)現(xiàn)液化后的砂土中剪應力和孔壓比之間有較好的線性關系。提出以下公式：

(8)

式中：τ為土體剪應力；pα為標準大氣壓；A和B為試驗參數(shù)。

式(8)可以反映土體發(fā)生液化時，強度和空隙水壓力的關系,但是不能反映變形和空隙水壓力之間的相互影響。

I. VARDOULAKIS等[32-33]選取孔隙率為耦合參數(shù)，建立了包含滲流侵蝕本構方程、混合滲流平衡方程和力平衡方程在內的數(shù)學模型，采用伽遼金有限元法模擬了油井出砂問題。J. BELL最早提出了滲流侵蝕本構方程：

(9)

式中：m為任意時刻土骨架相轉化為可動細顆粒相的速率；mcr為顆粒質量侵蝕速率；mdep為顆粒沉積速率；c為可動細顆粒濃度；qi為水相和可動顆粒相混合物的流量；φ為土體孔隙率；ρs為土骨架相密度。

管涌多相耦合多采用的有限差分和有限元模型或者不能體現(xiàn)管涌漸進性破壞特點或者不能考慮管涌發(fā)展過程中孔隙水壓力的變化，即不能考慮管涌發(fā)展過程的滲流侵蝕耦合效應。

滲流侵蝕的耦合機理將土體分為土骨架相、液相和液化顆粒相，基于多孔介質動力學理論建立三相質量守恒微分方程，再引入可以考慮土骨架的侵蝕導致液化顆粒相質量增加的滲流侵蝕本構方程，最后得出了一維三相滲流侵蝕耦合管涌數(shù)學模型。

L.SIBILLE等[34]使用離散元模擬土顆粒，格子玻爾茲曼方法模擬水流，然后在兩種方法之間建立耦合關系，從微觀角度研究了管涌現(xiàn)象中土顆粒與土骨架分離和土顆粒隨水流遷移兩個階段。分析時，土顆粒之間采用摩擦和黏結接觸本構，水流和土顆粒之間的作用，用水力剪力和來反映。

圖1 管涌產生機理Fig. 1 The mechanism of piping

圖2 管涌的多相耦合機理Fig. 2 The multiphase coupling mechanism of piping

注：其中u為位移；n為孔隙率；k為滲透系數(shù)；p為孔隙水壓力；q為滲透速度；c為細顆粒濃度；E為彈性模量；v為泊松比；φ為內摩擦角。

滲流侵蝕本構方程描述了液化細顆粒的流失量與滲流速度、管涌土體密度、孔隙率、液化細顆粒含量等物理量的關系，而現(xiàn)有管涌數(shù)學模型中常用的滲流侵蝕本構方程來源與石油工程中常用的砂巖的侵蝕方程，它與土木工程中常見的砂土，粉質土等的侵蝕性質不同，將其用于管涌數(shù)值模擬會產生誤差。

5 結論與展望

5.1 結 論

對土體流固耦合的4種分類方法，揭示每種理論的本質，有助于加深對流固耦合問題的認識，同時也指出了流固耦合存在的問題和發(fā)展方向。

1)計算域之間耦合的理論發(fā)展較早，計算簡單，但這是建立在過度簡化基礎上得出的，所以這種耦合精度較差。該種方法應進一步研究更復雜的耦合條件，例如不同介質的接觸條件。

2)基本未知數(shù)之間耦合的理論較為完善，耦合的計算模型與實際符合的較好，缺點是計算繁瑣，結合有限單元法可以在工程中得到充分的利用。

3)參數(shù)之間耦合能直觀的反應耦合現(xiàn)象，但是參數(shù)之間的耦合多是根據(jù)某一種土的實驗數(shù)據(jù)推導出來的，沒有普遍性，研究有待加深。

4)本構關系耦合是研究管涌、土體液化的理論基礎，但由于該理論的建立是在砂的基礎上，因此使用到粉土，黏性土中還需要進一步改進和修正。

5.2 展 望

從流固耦合現(xiàn)有理論中存在的問題中可以看出現(xiàn)有耦合理論還有很多不足，可做如下改進：

1)對計算域之間耦合的理論需要進一步改進，使流體域和固體域之間既滿足力的平衡條件也要滿足變形協(xié)調條件。

2)可以將參數(shù)之間的耦合應用到其他3種耦合形式中，例如Terzaghi一維固結理論只是滿足土體變形和液相排出體積的協(xié)調關系，而引入應力和滲透系數(shù)的關系后就可以考慮土體變形的非線性和滲流的非線性，使其更符合工程實際。

3)對本構關系耦合的研究較少，僅局限于對管涌、液化等問題的研究。以后，可以將本構關系的流固耦合推廣到普通的固結問題。