雨季夜巡的時(shí)候最怕看到擋墻泄水孔里流出渾水那種時(shí)候坡體內(nèi)部已經(jīng)在悄悄“鬧脾氣”了。真正讓我下定決心把非飽和入滲做進(jìn)數(shù)值模型里的是某次加固設(shè)計(jì)復(fù)核經(jīng)典條分法給出的安全系數(shù)接近1.18坡腳卻已經(jīng)出了錯(cuò)動(dòng)裂縫。反反復(fù)復(fù)核對(duì)參數(shù)之后問題并不在強(qiáng)度參數(shù)而在強(qiáng)度折減的前提——我把邊坡當(dāng)成滲流隨時(shí)達(dá)到飽和的土體完全沒有考慮雨水入滲過程中淺層負(fù)孔隙水壓力的快速耗散。這部分被短暫“借用”的強(qiáng)度一消失安全儲(chǔ)備立刻被透支。這篇文章就是圍繞 COMSOL 里“邊坡降雨不飽和條件下強(qiáng)度折減影響”的完整實(shí)現(xiàn)適合正在跟 COMSOL 死磕滲流-應(yīng)力耦合的同行也適合想把非飽和土力學(xué)落到具體算例的巖土工程師。你會(huì)看到從 Richards 方程、土水特征曲線參數(shù)、穩(wěn)態(tài)初場設(shè)置到強(qiáng)度折減的等效參數(shù)掃描以及最后如何判斷臨界折減系數(shù)。整個(gè)過程都是我實(shí)際跑通的一套流程不是教科書目錄。1. 勝負(fù)手在“吸力”非飽和強(qiáng)度到底被什么偷走1.1 負(fù)孔壓不是負(fù)擔(dān)而是一種隱形的黏聚力干燥邊坡在地下水位以上依然能立住很大程度靠的是非飽和區(qū)的負(fù)孔隙水壓力——也就是基質(zhì)吸力。土力學(xué)里常把吸力引起的抗剪強(qiáng)度增量寫成Δτ (u_a ? u_w) · tanφ?其中 u_a 是孔隙氣壓力u_w 是孔隙水壓力φ? 是吸力摩擦角。φ? 通常小于內(nèi)摩擦角 φ工程上很多取 φ 的一半甚至三分之一。你可以把它理解為吸力像往土里塞了無數(shù)根細(xì)小的“拉索”把顆粒相對(duì)位置額外固定了一層。降雨入滲時(shí)雨水把孔隙中的空氣逐步擠出負(fù)孔壓數(shù)值變小甚至變?yōu)檎龎豪饕桓傻?。這是強(qiáng)度折減背后真正的物理過程也正是“非飽和條件”和“飽和條件”最本質(zhì)的分水嶺。所以雨天才容易滑不是因?yàn)樗亓撕芏喽且驗(yàn)椤袄鳌痹诮饨?。這個(gè)解釋雖然粗糙但用來理解數(shù)值結(jié)果非常有效當(dāng)你看到安全系數(shù)驟降的時(shí)候先別急著算水壓力看看吸力場怎么退化。1.2 只算飽和水位線的模型錯(cuò)在哪傳統(tǒng)分析里很多做法把降雨影響等效為地下水位抬升然后據(jù)此計(jì)算孔隙水壓力增加。問題在于絕大多數(shù)淺層滑面發(fā)生在非飽和帶那里負(fù)孔壓的變化遠(yuǎn)比水位抬升劇烈。用飽和模型算往往出現(xiàn)兩種錯(cuò)一是整體安全系數(shù)偏高因?yàn)槟Ｐ蜎]有體現(xiàn)表層強(qiáng)度被削弱二是在坡體上部和坡頂處判不出塑性區(qū)而那恰是現(xiàn)場最容易出現(xiàn)張拉裂縫的位置。拿某均質(zhì)粉質(zhì)黏土邊坡算例來說同樣的強(qiáng)度參數(shù)飽和穩(wěn)態(tài)模型給 1.12而完整非飽和-瞬態(tài)模型在暴雨第 3 天給 0.94差別足以改變工程決策。我整理了一張對(duì)比表方便你直接感受差距對(duì)比項(xiàng)飽和穩(wěn)態(tài)近似非飽和瞬態(tài)入滲地下水位以上孔壓默認(rèn) 0不考慮吸力負(fù)孔壓按 SWCC 分布降雨影響路徑水位抬升吸力耗散 局部暫態(tài)飽和強(qiáng)度衰減來源正孔壓上升吸力項(xiàng)消失 正孔壓上升滑面位置預(yù)測偏深、偏坡腳淺層、坡腳至坡頂均可典型安全系數(shù)偏差偏高 10%~20%反映真實(shí)風(fēng)險(xiǎn)窗口這個(gè)代差就是本文所有數(shù)值設(shè)置的出發(fā)點(diǎn)。后面每一步操作目的都是把“非飽和”這三個(gè)字真實(shí)地放進(jìn)模型里。2. 讓雨水進(jìn)模型的正確姿勢從 Richards 方程到穩(wěn)態(tài)初場2.1 Van Genuchten 參數(shù)的物理解讀與取值COMSOL 里做非飽和滲流最省事的是用多孔介質(zhì)流物理接口里的 Richards 方程選項(xiàng)而不是手動(dòng)自定義 PDE因?yàn)樗詭?van Genuchten 土水特征曲線模型和相對(duì)滲透系數(shù)函數(shù)少出一堆手誤。Richards 方程簡化地寫成(C Se·S) ?h/?t ?·[?Ks·Kr(h)·?H] 0不用去背每一項(xiàng)你只需要知道土體儲(chǔ)水能力 C 和非飽和滲透系數(shù) Ks·Kr 的乘積決定了水分鋒面的推進(jìn)速度。Kr 隨吸力變化這是非飽和滲流和飽和滲流最大的區(qū)別——飽和模型里 Kr 恒等于 1非飽和模型里它可以在幾個(gè)數(shù)量級(jí)之間變動(dòng)。參數(shù)取值上以某山區(qū)公路邊坡的粉質(zhì)黏土為例我常用來作為基準(zhǔn)的一組值參數(shù)含義取值Ks飽和滲透系數(shù)5×10?? m/sθr殘余含水率0.03θs飽和含水率0.42α進(jìn)氣值相關(guān)參數(shù)1.6 1/mnSWCC 形態(tài)參數(shù)2.1m由 m1?1/n 得到0.52其中 α 影響進(jìn)氣值n 控制土水特征曲線的陡峭程度。如果手頭沒有實(shí)測建議參考同地區(qū)勘察報(bào)告不要把基坑經(jīng)驗(yàn)值直接套到邊坡。同一類型土壓實(shí)度和孔隙比不一樣SWCC 差別很大。2.2 穩(wěn)態(tài)初始孔壓場不跑這一步瞬態(tài)全都白搭很多同學(xué)急著拉時(shí)間軸一上來就給降雨邊界結(jié)果頭幾個(gè)時(shí)間步瘋狂不收斂。原因很簡單你給的初場和邊界條件自相矛盾。標(biāo)準(zhǔn)做法是先跑一個(gè)“無降雨穩(wěn)態(tài)”得到與實(shí)際地下水位對(duì)應(yīng)的初始孔壓分布。幾何模型我用的是坡高 10m、坡角約 43°、坡頂寬度 20m、坡腳延伸 15m 的均質(zhì)截面。地下水位設(shè)置在距坡腳底面約 8m 深處因此穩(wěn)態(tài)場里地下水面以下是正孔壓以上是負(fù)孔壓初始吸力隨高度增加。在 COMSOL 里操作時(shí)記得把 Richards 方程接口的研究步驟設(shè)為“穩(wěn)態(tài)”先算一遍然后在下一個(gè)瞬態(tài)研究中把穩(wěn)態(tài)解作為初值條件導(dǎo)入。這一步多花五分鐘能省下后面 50 次報(bào)錯(cuò)。2.3 降雨邊界的數(shù)量級(jí)陷阱mm/d 與 m/s 的換算降雨條件不是點(diǎn)擊“通量邊界”就完事單位換算首先要命。COMSOL 默認(rèn)用 SI 單位通量邊界寫 m/s氣象資料通常是 mm/d。換算關(guān)系1mm/d ≈ 1.157×10?? m/s暴雨 50mm/d 對(duì)應(yīng)約 5.8×10?? m/s強(qiáng)暴雨 200mm/d 對(duì)應(yīng)約 2.3×10?? m/s。如果氣象資料給 100mm/d把數(shù)字“100”直接填進(jìn)通量邊界相當(dāng)于每天下 1m 的雨邊坡不滑才怪。另一個(gè)容易踩的坑是降雨強(qiáng)度超過表層飽和滲透系數(shù)以后多余水量在模型里會(huì)積在表面制造假高壓。現(xiàn)實(shí)里這部分會(huì)形成坡面徑流數(shù)值上建議給入滲通量設(shè)上限比如q min(q_rain, Ks·β)β 取 1~2。這樣既不會(huì)讓模型被虛假積水撐爆又保留了降雨強(qiáng)度對(duì)入滲的控制作用。3. 強(qiáng)度折減在 COMSOL 里的非標(biāo)準(zhǔn)落地參數(shù)掃描代替“一鍵折減”3.1 為什么 COMSOL 沒有“自動(dòng)折減”按鈕用慣了有限差分強(qiáng)度折減程序的人會(huì)不習(xí)慣那邊一個(gè)命令下去自動(dòng)二分折減系數(shù)COMSOL 里并沒有這么個(gè)按鈕。沒說不能做只是需要把它拆成參數(shù)化掃描。思路是把材料參數(shù)里的有效黏聚力 c′ 和有效內(nèi)摩擦角 φ′ 定義成全局參數(shù) FS 的表達(dá)式然后讓求解器對(duì) FS 從 1.0 掃描到 2.0。具體做法是全局參數(shù)里寫 FS1材料參數(shù)里把 c′ 替換成 c0/FS把 tanφ′ 替換成 tanφ0/FS再加一個(gè)輔助掃描。FS 從 1.0 起步每步遞增 0.1接近臨界步長加密到 0.02一共算十幾步足夠捕捉突變。掃描時(shí)建議開啟“繼續(xù)”選項(xiàng)讓每個(gè) FS 的求解都從上一個(gè) FS 的解出發(fā)。這樣做的好處是塑性發(fā)展歷史連續(xù)收斂速度明顯加快物理上也更合理——真實(shí)工程的荷載增加本來就是連續(xù)的。3.2 折減公式的非飽和變形吸力項(xiàng)到底折不折這是和純飽和模型最大的不同點(diǎn)。非飽和強(qiáng)度如果寫成τ_f c″/FS (σ_n ? u_a) · tanφ′/FS (u_a ? u_w) · tanφ?/FS那就是包絡(luò)式寫法把吸力項(xiàng)也放進(jìn)了折減。另一種做法是把吸力項(xiàng)當(dāng)作外環(huán)境效應(yīng)不參與折減折減只作用于原位強(qiáng)度τ_f (c″ (σ_n ? u_a) · tanφ′)/FS (u_a ? u_w) · tanφ?哪種對(duì)取決于研究目的。影響研究通常兩種都跑工程復(fù)核我會(huì)偏向第一種更保守、也更容易被評(píng)審接受。實(shí)際算例中當(dāng) φ? 取 φ 的三分之一時(shí)兩種方案的安全系數(shù)差異大約在 5%~8%趨勢一致但數(shù)值略有區(qū)別。如果你在報(bào)告里寫“吸力項(xiàng)不參與折減”一定要把這一條單列出來別混在所有結(jié)果里否則別人復(fù)現(xiàn)數(shù)據(jù)時(shí)會(huì)對(duì)不上。3.3 從塑性區(qū)貫通到位移拐點(diǎn)臨界折減判據(jù)怎么定掃描完成后怎么判斷“臨界折減系數(shù)”標(biāo)準(zhǔn)判據(jù)有兩個(gè)一是塑性應(yīng)變等值線從坡腳向坡頂貫通形成連續(xù)的剪切滑移通道取剛好貫通的 FS 近似臨界值。這個(gè)判據(jù)直觀但需要人為判斷“貫通”的瞬刻有時(shí)候塑性區(qū)看起來接近貫通實(shí)際差一口氣。二是看坡頂位移與 FS 曲線位移突然從線性段進(jìn)入陡增段的那個(gè) FS就是臨界折減系數(shù)。位移拐點(diǎn)更客觀適合批量處理。實(shí)操上我建議兩個(gè)判據(jù)一起用交叉驗(yàn)證。COMSOL 里畫出“坡頂某點(diǎn)的位移-FS”曲線不用額外編程后處理里取點(diǎn)就行。判斷的時(shí)候我習(xí)慣把塑性區(qū)云圖疊加在位移曲線上如果拐點(diǎn) FS 對(duì)應(yīng)塑性區(qū)正好剛貫通說明結(jié)果非常“干凈”如果兩者對(duì)不上多半是網(wǎng)格不夠密或者掃描步長太大需要回去加密重算。4. 算例說話一個(gè)均質(zhì)粉質(zhì)黏土邊坡降雨 7 天的安全系數(shù)演化4.1 基準(zhǔn)工況FS 隨降雨歷時(shí)的變化把前面說的一整套參數(shù)放到算例里。基準(zhǔn)工況Ks5×10?? m/sα1.6 1/m初始水位 8m降雨強(qiáng)度 80mm/d連續(xù)下 3 天、停 4 天。計(jì)算結(jié)果如下時(shí)間節(jié)點(diǎn)臨界折減系數(shù) FS降雨前1.18降雨 12h 后1.07降雨 24h 后1.01降雨 72h 后0.94停雨后第 2 天0.97停雨后第 4 天1.02降雨結(jié)束之后FS 從 0.94 只恢復(fù)到 1.02遠(yuǎn)沒有回到初始的 1.18。這說明降雨影響不只是“下的時(shí)候危險(xiǎn)”停雨后的滯后效應(yīng)同樣要命。如果只看降雨結(jié)束時(shí)刻會(huì)低估邊坡實(shí)際的風(fēng)險(xiǎn)窗口。物理解釋也不復(fù)雜入滲前鋒掃過淺層吸力銳減但孔壓恢復(fù)依賴土體排水和蒸發(fā)滲透系數(shù)越低恢復(fù)越慢同時(shí)坡腳處因?yàn)閼?yīng)力集中最先出現(xiàn)塑性應(yīng)變即使整體沒有貫通位移已經(jīng)不可逆。4.2 滲透系數(shù)、吸力參數(shù)與初始水位的影響規(guī)律只跑一條曲線不夠影響研究需要做參數(shù)敏感性。我列了一個(gè) 5 組工況的矩陣每組在相同降雨條件下跑完 7 天工況Ks (m/s)α (1/m)初始水位 (m)3d FS7d FS基準(zhǔn)5×10??1.680.941.02高滲透5×10??1.681.081.14低滲透5×10??1.680.880.91緩吸力5×10??0.681.031.10淺水位5×10??1.630.730.78幾條規(guī)律值得記住第一初始水位淺的邊坡風(fēng)險(xiǎn)最大安全系數(shù)直接從安全區(qū)進(jìn)入危險(xiǎn)區(qū)而且 7 天恢復(fù)量很小。這類邊坡雨季前就應(yīng)該重點(diǎn)關(guān)注。第二低滲透性土反而更危險(xiǎn)表面容易形成暫態(tài)飽和帶降水進(jìn)得去出不來峰值強(qiáng)度損失大。第三α 越小SWCC 越緩吸力維持能力越強(qiáng)FS 下降越慢。這套規(guī)律提醒做影響研究的人敏感性分析至少要做 2×2 矩陣別只改一個(gè)參數(shù)單點(diǎn)對(duì)比。二維邊坡模型跑一次不算貴但結(jié)論的說服力完全不同。4.3 與條分法的交叉驗(yàn)證數(shù)值模型最怕自說自話。把相同強(qiáng)度參數(shù)、相同孔壓場導(dǎo)出到經(jīng)典簡化 Bishop 條分程序復(fù)核關(guān)鍵工況的安全系數(shù)差在 0.02~0.05 以內(nèi)趨勢完全一致。差值的來源主要是條分法對(duì)縱向應(yīng)力積分做了簡化以及二維平面應(yīng)變假設(shè)。這個(gè)操作能顯著提高結(jié)果說服力特別是拿到評(píng)審面前。我的習(xí)慣是先把非飽和數(shù)值模型的 FS 曲線畫出來再把條分法算的幾個(gè)狀態(tài)點(diǎn)疊到同一張圖上。如果兩者只在 0.03 以內(nèi)錯(cuò)動(dòng)那這組結(jié)果基本可信如果差到 0.1 以上優(yōu)先排查滲流場有沒有算對(duì)而不是力學(xué)參數(shù)。5. 收斂、網(wǎng)格和判據(jù)我只留下來三條真管用的經(jīng)驗(yàn)5.1 不收斂時(shí)先自問三件事第一件事初場穩(wěn)態(tài)算過沒有。瞬態(tài)分析直接帶降雨邊界10 次里有 8 次前幾個(gè)時(shí)間步瘋狂振蕩。第二件事時(shí)間步長是不是給太兇。Richards 方程是很頑固的非線性項(xiàng)尤其是 SWCC 斜率大的區(qū)域建議求解時(shí)間序列初期用對(duì)數(shù)分布從 1e-2 小時(shí)到 12 小時(shí)而不是從 0 直接跨到 24 小時(shí)。第三件事塑性更新和滲流更新是否同步。如果滲流和力學(xué)在同一研究里做全耦合計(jì)算量爆炸且收斂極差。我更推薦先算滲流把孔壓場導(dǎo)出成“解”再在力學(xué)研究里按時(shí)間點(diǎn)插入。兩步解耦在我們這套非飽和強(qiáng)度折減流程里效率提高了一個(gè)數(shù)量級(jí)。5.2 邊界層網(wǎng)格與時(shí)間步長組合網(wǎng)格不是越密越好而是“該密的地方密到可控”。坡面入滲從邊界開始表層 1m 范圍內(nèi)的孔壓梯度最陡我習(xí)慣在坡表加邊界層網(wǎng)格首層厚度 0.1m、偏置系數(shù) 1.2總網(wǎng)格控制在 1 萬到 3 萬之間坡體內(nèi)部用 1m 的三角形單元坡腳處再局部加密到 0.2m。以基準(zhǔn)工況為例滲流瞬態(tài)大約算半小時(shí)加上 10 次參數(shù)掃描的力學(xué)求解大概兩個(gè)多小時(shí)普通工作站就能跑完。時(shí)間步長配合網(wǎng)格有一條鐵律表層單元越薄初期時(shí)間步就要越小否則高頻振蕩會(huì)在孔壓云圖上留下棋盤狀跳躍后處理根本沒法看。5.3 判斷結(jié)果可信的最小自檢清單最后梳理一個(gè)我每次算完都會(huì)過的清單你也可以直接復(fù)制當(dāng)模板初始穩(wěn)態(tài)孔壓場的地下水位線是否與設(shè)定一致坡頂負(fù)壓數(shù)值是否在合理范圍降雨期間表層負(fù)孔壓是否隨時(shí)間單調(diào)下降有沒有前幾步回彈的怪相塑性區(qū)是否從坡腳附近先出現(xiàn)再逐步向坡頂延伸如果先從坡頂出現(xiàn)說明網(wǎng)格或邊界有原則性錯(cuò)誤臨界折減系數(shù)與條分法的差距是否在 0.05 以內(nèi)報(bào)告里所有 FS 值是否注明了對(duì)應(yīng)降雨時(shí)間和折減策略吸力項(xiàng)是否參與折減。通常前三條有問題先回模型后兩條有問題先回工程判斷。這五條逐個(gè)過完我才敢把結(jié)果放進(jìn)正式計(jì)算書。最后補(bǔ)一句個(gè)人體會(huì)吧。如果只能記住一件事我建議做這個(gè)方向的人把 φ? 的取值和吸力項(xiàng)折減策略提前寫進(jìn)報(bào)告因?yàn)樗苯記Q定 FS 的數(shù)值口徑。另一個(gè)人可以馬上用的習(xí)慣是把幾何、邊界層網(wǎng)格、Richards 接口和參數(shù)掃描模板存成模型模板下次換土性參數(shù)改三處數(shù)值就能跑出新方案。我近幾次做降雨邊坡影響評(píng)估都靠這個(gè)模板省掉的重復(fù)勞動(dòng)遠(yuǎn)大于初期建模板的時(shí)間這算是能給同行最實(shí)際的一條建議。