土石壩滲流的流固耦合方法計算分析

代凌輝，張 玓，楊衛(wèi)元

（1.黃河水利職業(yè)技術(shù)學院，河南 開封 475004；2.鶴慶縣水利水電勘測設(shè)計隊，云南 鶴慶 671500）

0 引言

在土石壩的設(shè)計中，滲流問題對壩體的穩(wěn)定性至關(guān)重要，也極大地影響著構(gòu)筑物的安全和造價。 在進行排水和防滲設(shè)施設(shè)計或進行壩體穩(wěn)定性計算時，一般都需要先確定出滲流的自由面[1～2]。 在目前的分析方法中，由于有限元方法能夠有效地處理復(fù)雜的邊界條件、材料的非均勻性、材料的各向異性，并能方便地求解三維問題，在工程設(shè)計中被廣泛使用。 傳統(tǒng)的滲流有限元求解方法的難點是自由面（浸潤面）的位置確定，因為滲流區(qū)域是未知的，需要通過迭代求解。 使用變網(wǎng)格迭代法，在每次迭代計算時，都要確定自由面的位置，并根據(jù)自由面的位置調(diào)整滲流網(wǎng)格。 這樣，在計算過程中，不但需要重新形成和分解總體傳導矩陣，耗費大量的計算時間，而且在自由面附近的單元也可能出現(xiàn)畸形，使得求解失真[3～4]。 目前，許多研究者采用固定網(wǎng)格的方法進行研究。 該方法按照非飽和土力學理論，將整個區(qū)域作為分析區(qū)域，浸潤面取孔隙水壓力為零的地方，為問題的求解提供了很大的便利，并保證了較高的精度[5]。 本文試采用基于固定網(wǎng)格的有限元流固耦合法進行分析。

流固耦合是研究變形固體在流場作用下的各種行為，以及固體位行對流場影響的一種計算方法，其重要特征是兩相介質(zhì)之間的相互作用。 變形固體在流體載荷作用下會產(chǎn)生變形或運動，變形或運動又反過來影響滲流，從而改變流體載荷的分布和大小。

1 土石壩非飽和滲流的典型邊界條件

土石壩非飽和滲流在進行流固耦合計算時的典型邊界條件如圖1 所示。 其各邊的邊界條件為[5]：（1）S1是上游壩面水位以下部分，邊界條件為上游已知總水頭（此處的總水頭等于H1）。在進行有限元流固耦合分析時，只需給出邊界上的孔隙水壓力即可。 （2）S2是下游壩面水位以下部分，邊界條件為下游已知總水頭（此處的總水頭等于H2）。 同樣，在進行有限元流固耦合分析時，只需要給出邊界上的孔隙水壓力。 （3）S3是壩基不透水基巖，為不透水邊界條件，通過該邊界的流量為零。 （4）S4為浸潤線，其位置是未知的。 該邊界上孔隙水壓力為零，即總水頭等于位置水頭z。 在傳統(tǒng)的滲流分析中，將該面也看成不透水邊界，即滲流只在飽和區(qū)內(nèi)發(fā)生，浸潤面是最上方的流線。 但在非飽和滲流分析中，無須如此設(shè)置。 （5）S5為自由滲出段邊界，孔隙水壓力為零，且滲流只能沿著下游坡面進行。 在有限元軟件中，該種邊界條件只允許孔隙水從分析區(qū)域中流出，而不允許水流進入。 即假設(shè)當邊界面上的孔隙水壓力為正時，孔隙流體的流速與孔隙水壓力成正比；當孔隙水壓力為負值（吸力）時，流速限定為零，即流體不會進入到分析區(qū)域內(nèi)部。 具體表達式為：

式中：vn為滲流流速，ks為滲流系數(shù)，uw為孔隙水壓力。

圖1 土石壩非飽和滲流典型邊界條件Fig.1 Unsaturated seepage typical boundary conditions of earth dam

2 土石壩流固耦合計算

2.1 均質(zhì)土石壩流固耦合計算

某均質(zhì)土石壩壩高為15 m，壩頂寬度為4 m，上、下游坡比均為1∶2，壩底長度為64 m，上游水位高出建基面13 m，下游水位與建基面平齊。 壩殼料土體的滲透系數(shù)為1.0×10-7m/s。 建立有限元模型，共劃分306 個節(jié)點、272 個單元。 按照非飽和滲流典型邊界條件，賦予有限元模型相應(yīng)的邊界條件。 計算完成后，繪制孔隙水壓力等值線圖（如圖2 所示）。 從計算結(jié)果中提取迎水面上各節(jié)點的滲透流量值，相加之后得到該計算模型的滲透流量，為2.55×10-7m3/(s·m)。 該值與水力學方法得到的結(jié)果相比，相對誤差小于2%，由此證明流固耦合方法的計算結(jié)果比較準確。 從孔隙水壓力等值線圖中可以得出，孔隙水壓力為0 的線即為浸潤線，滲流逸出點的位置位于下游壩坡6 m 高處（如圖2 所示）。

圖2 均質(zhì)壩算例孔隙水壓力等值線圖（單位：kPa）Fig.2 Pore water pressures isoline of the homogeneous dam (Unit: kPa)

2.2 黏土心墻土石壩流固耦合計算

某黏土心墻土石壩壩高為127.5 m，壩頂寬為11 m，上游邊坡系數(shù)為1∶2，下游邊坡系數(shù)為1∶1.75，上游圍堰高為40 m，下游排水棱體高為30 m。 以建基面為相對0 高程面，上游水位117.5 m，下游水位10 m。 黏土心墻頂高程與壩頂高程一致，心墻頂寬11 m，底寬87.5 m。壩體的斷面如圖3 所示。壩體各分區(qū)的滲透系數(shù)見表1。 采用有限元流固耦合方法進行計算和分析，建立有限元模型，共劃分1 062 個節(jié)點，975 個單元。 按照非飽和滲流典型邊界條件，賦予有限元模型相應(yīng)的邊界條件。

圖3 黏土心墻壩橫斷面圖Fig.3 Cross-section of the clay core wall dam

表1 黏土心墻壩算例各分區(qū)的滲透系數(shù)Table 1 Permeability coefficients of each partition of the clay core wall dam

計算完成后，繪制該黏土心墻壩的孔隙水壓力等值線圖（如圖4 所示）。 從計算結(jié)果中提取迎水面上各節(jié)點的滲透流量值，相加之后，得到該計算模型的滲透流量，為5.19×10-6m3/(s·m)。從孔隙水壓力等值線圖可以得出，孔隙水壓力為0 的線為浸潤線?？紫端畨毫Φ戎稻€在黏土心墻部位急劇下降，說明該壩的黏土心墻防滲效果明顯。由此結(jié)果可以得出，對于壩體有不同分區(qū)、滲透系數(shù)不同的非均質(zhì)壩，有限元流固耦合方法依然適用，而且便捷、高效、結(jié)果準確。

圖4 黏土心墻壩孔隙水壓力等值線圖（單位：kPa）Fig.4 Pore water pressures isoline of the clay core wall dam (Unit: kPa)

3 結(jié)語

通過對兩個算例進行計算，主要得出以下結(jié)論：（1） 賦予有限元模型正確的邊界條件是得出準確結(jié)果的前提。（2）把計算結(jié)果中土石壩迎水面上各節(jié)點的滲透流量值相加，可得出土石壩整體的滲透流量。（3）根據(jù)計算結(jié)果，可以繪制出土石壩的孔隙水壓力等值線圖。 圖中，孔隙水壓力為0 的線即為浸潤線。對于下游逸出點高于下游水位的情況，還可由此確定逸出點的位置。 （4）對于壩體有不同分區(qū)、滲透系數(shù)不同的非均質(zhì)壩（如黏土心墻壩），有限元流固耦合方法依然適用，而且便捷、高效、結(jié)果較準確。非均質(zhì)壩的孔隙水壓力等值線會在心墻部位急劇下降，由此為檢驗心墻設(shè)計方案的合理性提供了方法和依據(jù)。

本文結(jié)果證明了采用有限元流固耦合方法計算土石壩非飽和滲流問題的可行性，可為相關(guān)類似工程的分析提供參考和借鑒。

[1] 金生，耿艷芬，王志力. 利用飽和-非飽和滲流模型計算壩體自由面滲流 [J]. 大連理工大學學報，2004，44(1):110-113.

[2] 牛文杰，葉為民. 不透水地基上均質(zhì)土壩滲流自由面的有限元計算[J] .巖土力學，2007（增刊）:375-378.

[3] 徐千軍，張建紅. 確定穩(wěn)定滲流自由面位置的一種簡便方法[J]. 水動力學研究與進展，1999, 14(4):418-423.

[4] 殷宗澤. 土工原理[M]. 北京:中國水利水電出版社，2007：171-189.

[5] 費康，張建偉. ABAQUS 在巖土工程中的應(yīng)用[M]. 北京：中國水利水電出版社，2010：205-206.