干數(shù)值模擬這行當(dāng)最讓人頭禿的往往不是力學(xué)推導(dǎo)本身而是建模那一關(guān)。常規(guī)有限元里裂縫面必須跟單元邊界重合網(wǎng)格拓?fù)涞眠w就裂縫幾何。遇到一條彎彎曲曲的裂縫倒還好頂多多切幾刀可真到了多裂縫交叉、分叉、隨機(jī)分布的場景網(wǎng)格重構(gòu)的代價(jià)能把人逼瘋。擴(kuò)展有限元XFEM就是沖著這個(gè)痛點(diǎn)來的它還允許裂縫“穿”過單元內(nèi)部位移不連續(xù)通過額外的富集項(xiàng)表達(dá)網(wǎng)格不用再貼著裂縫走。這個(gè)思路特別對搞水力壓裂、石油工程、地質(zhì)力學(xué)數(shù)值模擬的朋友胃口因?yàn)樗烊贿m合處理“裂縫網(wǎng)絡(luò)裂縫滲透率”這種聯(lián)動(dòng)問題。這篇文章我會(huì)把XFEM從原理到C實(shí)現(xiàn)的關(guān)鍵環(huán)節(jié)拆開講尤其是裂縫與滲透率聯(lián)動(dòng)的處理給出能落地的代碼片段和實(shí)操建議適合想上手XFEM的入門者也適合已經(jīng)在用但被各種玄學(xué)bug困擾的同行。1. 為什么裂縫模擬繞不開擴(kuò)展有限元1.1 常規(guī)有限元的“網(wǎng)格牢籠”傳統(tǒng)有限元處理裂縫最粗暴的辦法是把裂縫面當(dāng)成兩個(gè)互相接觸的邊界網(wǎng)格必須在裂縫處斷開讓兩側(cè)節(jié)點(diǎn)獨(dú)立。這個(gè)做法原理上沒問題但實(shí)操極其痛苦。一條貫穿多個(gè)單元的裂縫你要在幾何建模階段就把網(wǎng)格“切”得跟裂縫完全一致等于讓裂縫牽著網(wǎng)格走。更麻煩的是一旦裂縫要擴(kuò)展整個(gè)局部網(wǎng)格就得重構(gòu)舊網(wǎng)格的場變量還要插值映射到新網(wǎng)格上精度損失和實(shí)現(xiàn)復(fù)雜度同時(shí)爆炸。我在實(shí)際項(xiàng)目里試過用界面單元cohesive element做裂縫模擬它比純斷開網(wǎng)格溫和一些裂縫路徑可以預(yù)先設(shè)置但也只適合預(yù)設(shè)裂縫面或者裂縫沿著單元邊界走的情形。一旦裂縫路徑是材料響應(yīng)的一部分比如壓裂過程中裂縫隨機(jī)轉(zhuǎn)向網(wǎng)格就得隨時(shí)待命計(jì)算成本大部分耗在重構(gòu)上而不是求解上。這時(shí)候你會(huì)特別渴望一種方法讓裂縫從網(wǎng)格的“緊身衣”里解脫出來。1.2 擴(kuò)展有限元的破局思路XFEM的核心思想說白了就一句話位移場里額外加幾項(xiàng)“補(bǔ)丁函數(shù)”讓不連續(xù)天然地存在于單元內(nèi)部。網(wǎng)格照常劃分完全無視裂縫位置裂縫在哪由水平集函數(shù)定位裂縫帶來的位移跳躍用Heaviside富集函數(shù)壓在裂縫穿過的那些節(jié)點(diǎn)上裂尖的應(yīng)力奇異性用一組分支函數(shù)去逼近。這聽起來像是偷懶其實(shí)是把幾何復(fù)雜度換成了自由度復(fù)雜度。力學(xué)上位移逼近多出來幾個(gè)未知量編程上自由度編號(hào)、單元積分要單獨(dú)照顧富集區(qū)域。但整體而言網(wǎng)格生成和裂縫建模解耦了你可以在規(guī)則網(wǎng)格上擺任意形狀的裂縫這在工程上帶來的便利是巨大的。對我這種經(jīng)常要試幾十種裂縫分布方案的人來說XFEM把“改裂縫形貌”從“重新畫網(wǎng)格”降級(jí)成了“改幾個(gè)水平集數(shù)值”效率提升是數(shù)量級(jí)的。1.3 裂縫與滲透率的聯(lián)動(dòng)問題把裂縫和滲透率放在一起是這類項(xiàng)目最值的部分。巖石基質(zhì)本身滲透率很低流體主要是沿著裂縫網(wǎng)絡(luò)流動(dòng)的。裂縫開度哪怕只有零點(diǎn)幾毫米導(dǎo)流能力也比基質(zhì)高出好幾個(gè)數(shù)量級(jí)。水力壓裂的本質(zhì)就是人為制造高導(dǎo)流裂縫把儲(chǔ)層里的油氣解放出來。要算準(zhǔn)這個(gè)“解放”過程至少需要兩層信息第一層是應(yīng)力場打出來的裂縫幾何狀態(tài)包括裂縫位置、開度分布第二層是裂縫幾何狀態(tài)到滲透率張量的映射關(guān)系。這兩層是耦合的應(yīng)力加載改變開度開度改變滲透率滲透率改變孔隙壓力孔隙壓力反過來又改變有效應(yīng)力。所以一個(gè)真正能用的代碼不能只算位移場還要能把位移結(jié)果“翻譯”成滲透率參數(shù)再接一個(gè)流動(dòng)求解器。這篇文章后面會(huì)逐步把這套流程打通。2. 核心原理拆解把“不連續(xù)”裝進(jìn)單元2.1 Heaviside富集一條裂縫的兩張臉從數(shù)學(xué)上講裂縫給位移場帶來的最直接變化就是“跳變”裂縫面上某一點(diǎn)的位移兩側(cè)可以不一樣。常規(guī)有限元的位移逼近是連續(xù)的它在單元內(nèi)部不可能表達(dá)跳變。XFEM的做法是在裂縫穿過的那些節(jié)點(diǎn)上除了常規(guī)形函數(shù)對應(yīng)的自由度再加一個(gè)帶Heaviside函數(shù)的富集項(xiàng)。整體位移逼近寫成這個(gè)形式u_h(x) Σ N_i(x) u_i Σ N_j(x) H(φ(x)) a_j ...第一項(xiàng)是標(biāo)準(zhǔn)有限元部分第二項(xiàng)是Heaviside富集。這里的H(φ)其實(shí)是個(gè)符號(hào)函數(shù)φ 0的裂縫一側(cè)取1另一側(cè)取-1。實(shí)際編程的時(shí)候我們把富集項(xiàng)寫成 N_j(x) (H(φ) - H(φ_j))減掉形函數(shù)在節(jié)點(diǎn)處的值是為了保證富集自由度a_j不影響節(jié)點(diǎn)處的位移值這樣邊界條件處理起來就跟標(biāo)準(zhǔn)有限元一樣清爽。這里有個(gè)容易被忽略的細(xì)節(jié)當(dāng)某個(gè)富集節(jié)點(diǎn)所在的單元里裂縫的幾何信息發(fā)生變化比如裂縫擴(kuò)展了一截節(jié)點(diǎn)是否富集、富集值是否要更新都得跟著變。這就是動(dòng)態(tài)富集也是程序中最容易出bug的地方之一。我最初實(shí)現(xiàn)的時(shí)候是在每次增量步開頭重新掃描一遍所有單元與裂縫的位置關(guān)系代價(jià)稍高但邏輯清晰實(shí)測穩(wěn)定。2.2 裂尖分支函數(shù)奇異性怎么表達(dá)裂縫尖端是另一個(gè)麻煩。彈性力學(xué)解析解告訴我們裂尖附近應(yīng)力場有r^{-1/2}量級(jí)的奇異性位移場在裂尖呈√r的漸近行為。如果你只用Heaviside富集去描述一條裂縫裂尖那個(gè)單元還是處理不了裂縫面在那里“剎車”了但應(yīng)力奇異性沒人管。解決方法是給裂尖所在的節(jié)點(diǎn)再疊加一組分支富集函數(shù)經(jīng)典的四個(gè)分支函數(shù)長這樣{ √r·sin(θ/2), √r·cos(θ/2), √r·sin(θ/2)·sinθ, √r·cos(θ/2)·sinθ }其中r和θ是裂尖局部極坐標(biāo)。實(shí)際實(shí)現(xiàn)時(shí)注意θ的原子取值范圍是(-π, π]asin、atan2函數(shù)的返回區(qū)間要處理對否則富集函數(shù)會(huì)在裂縫面上出現(xiàn)不期望的間斷。使用裂尖富集是XEFM實(shí)現(xiàn)中進(jìn)階玩法。如果你的重點(diǎn)是裂縫網(wǎng)絡(luò)滲透率分析裂縫尖端區(qū)域相對整個(gè)模型占比不大裂尖富集的精度收益主要體現(xiàn)在應(yīng)力強(qiáng)度因子K_I、K_II的提取上。如果只關(guān)心裂縫開度和滲透率分布穩(wěn)定捕捉位移跳躍比追求裂尖奇異性精度更重要——這點(diǎn)我推薦初學(xué)階段先把Heaviside富集做到萬無一失再加裂尖富集。2.3 水平集函數(shù)幾何信息與富集區(qū)判定剛才代碼里出現(xiàn)的φ就是水平集函數(shù)。它在XFEM里的角色是“裂縫幾何的數(shù)學(xué)翻譯官”。對一條裂縫我們定義符號(hào)距離函數(shù)φ(x)裂縫面上φ0裂縫一側(cè)φ0另一側(cè)φ0。對每個(gè)節(jié)點(diǎn)存下φ值就能判斷裂縫在哪個(gè)單元內(nèi)部穿過以及節(jié)點(diǎn)位于裂縫哪一側(cè)。但只有φ還不夠定位裂尖。還要一個(gè)裂縫面切向的水平集函數(shù)ψψ0的等值線過裂尖且垂直于裂縫面。有了φ和ψ裂尖位置就是兩個(gè)水平集同時(shí)為0的點(diǎn)。富集節(jié)點(diǎn)區(qū)域的判定也要靠它們某個(gè)單元內(nèi)有節(jié)點(diǎn)φ符號(hào)相反說明裂縫穿過該單元節(jié)點(diǎn)做Heaviside富集某個(gè)節(jié)點(diǎn)滿足 r √(φ2ψ2) ≤ r_enrich就能納入裂尖富集范圍。這個(gè)r_enrich一般取裂尖周圍一兩層單元的尺寸太大富集節(jié)點(diǎn)過多矩陣變得病態(tài)太小分支函數(shù)逼不出奇異應(yīng)力場。水平集函數(shù)還有一個(gè)好處支持復(fù)雜裂縫網(wǎng)絡(luò)。多條裂縫就存多組水平集裂縫分叉、相交也不過是多個(gè)水平集函數(shù)的布爾組合問題。幾何操作從“重構(gòu)網(wǎng)格”變成“更新標(biāo)量場”這是XFEM能應(yīng)付復(fù)雜裂縫的根本原因。3. C實(shí)現(xiàn)架構(gòu)從零搭建一個(gè)XFEM求解器3.1 數(shù)據(jù)結(jié)構(gòu)節(jié)點(diǎn)、單元、裂縫對象寫XFEM程序最忌諱的就是沿用標(biāo)準(zhǔn)有限元的“節(jié)點(diǎn)—單元”兩件套往里面硬塞富集信息。我建議一開始就把數(shù)據(jù)結(jié)構(gòu)設(shè)計(jì)成三件套節(jié)點(diǎn)、單元、裂縫。每個(gè)節(jié)點(diǎn)除了坐標(biāo)必須帶自己的水平集值每條裂縫一個(gè)φ值和富集標(biāo)記每個(gè)單元要記錄它可能包含的裂縫段和裂尖裂縫對象獨(dú)立管理幾何和開度歷史。一個(gè)最簡數(shù)據(jù)結(jié)構(gòu)大概長這樣struct Node { int id; double x, y; bool he_rich; // 是否被Heaviside富集 bool tip_rich; // 是否被裂尖富集 std::vectordouble phi; // 與每條裂縫對應(yīng)的水平集值 double psi; // 裂尖水平集按裂縫編號(hào)存多個(gè)也行 double ux, uy; // 位移解 }; struct CrackSegment { int id; std::vectorstd::arraydouble,2 pts; // 折線段 std::vectordouble opening; // 各分段開度 }; struct Element { int id; std::vectorint nodes; bool crossed; // 有裂縫穿過 bool has_tip; // 含裂尖 std::vectorint crack_ids; // 與哪些裂縫相交 };這種結(jié)構(gòu)的好處是富集信息跟著節(jié)點(diǎn)走裂縫幾何單獨(dú)存在單元只做交集判定。后處理算開度、算滲透率時(shí)直接從CrackSegment里取數(shù)據(jù)就行不需要每算一次都重新遍歷網(wǎng)格。3.2 自由度裝配富集項(xiàng)如何編號(hào)XFEM程序里最勸退新手的環(huán)節(jié)恐怕就是自由度裝配。常規(guī)有限元一個(gè)節(jié)點(diǎn)兩個(gè)位移自由度二維問題很好編號(hào)。XFEM的節(jié)點(diǎn)如果被富集了每多一個(gè)富集項(xiàng)就多一組自由度。我的裝配策略是這樣的先統(tǒng)計(jì)每個(gè)節(jié)點(diǎn)的富集類型Heaviside富集追加一個(gè)自由度對裂尖富集追加四個(gè)自由度對。然后遍歷節(jié)點(diǎn)給所有自由度排全局編號(hào)。編號(hào)順序本身沒講究但一定要保證“每個(gè)節(jié)點(diǎn)的所有自由度編號(hào)連續(xù)”這樣后面組裝單元?jiǎng)偠染仃嚂r(shí)取局部到全局的映射表才順手。int globalDof 0; for (auto nd : nodes) { nd.dof_std_x globalDof; // 標(biāo)準(zhǔn)x自由度 nd.dof_std_y globalDof; if (nd.he_rich) { nd.dof_he_x globalDof; nd.dof_he_y globalDof; } if (nd.tip_rich) { for (int a 0; a 4; a) { nd.dof_tip_x[a] globalDof; nd.dof_tip_y[a] globalDof; } } }這一步做對了后面組裝剛度矩陣就只是“對號(hào)入座”。很多所謂“矩陣奇異、無異常退出”的bug根源都是裝配時(shí)自由度編號(hào)重疊或者遺漏排查起來讓人懷疑人生。裝配完可以寫個(gè)小函數(shù)打印總的自由度數(shù)量跟手算的對一遍能在早期攔截大量問題。3.3 單元?jiǎng)偠染仃嚺c整體求解組裝單元?jiǎng)偠染仃嚨牧鞒谈鷺?biāo)準(zhǔn)有限元一樣先算形函數(shù)導(dǎo)數(shù)、雅可比再疊加重積。唯一不同的是被裂縫穿過的單元形函數(shù)空間里多了富集項(xiàng)的導(dǎo)數(shù)裂尖單元?jiǎng)t要算裂尖局部坐標(biāo)由水平集值轉(zhuǎn)換而來再求分支函數(shù)的導(dǎo)數(shù)。以二維線彈性問題為例B矩陣要擴(kuò)展成包含標(biāo)準(zhǔn)項(xiàng)和富集項(xiàng)的組合。對被裂縫穿過的單元B矩陣的額外列來自N_j乘以H在全局坐標(biāo)下的導(dǎo)數(shù)。由于H在裂縫面兩側(cè)是常數(shù)±1導(dǎo)數(shù)其實(shí)就是零——這就是所謂“Heaviside項(xiàng)的應(yīng)變貢獻(xiàn)只在裂尖區(qū)域有非零值”的由來但在富集自由度對應(yīng)節(jié)點(diǎn)的整體方程里形函數(shù)導(dǎo)數(shù)仍然要按N_j本身來算不能因?yàn)镠導(dǎo)數(shù)為零就整個(gè)消掉這一點(diǎn)困擾我好久。說白了H的作用是把富集自由度的“作動(dòng)量”集中在裂縫兩側(cè)而不是在裂縫內(nèi)部制造應(yīng)變。線性方程組的求解我強(qiáng)烈建議直接用成熟的稀疏直接求解器別再自己造輪子。對二維、幾千個(gè)自由度以內(nèi)的模型用Eigen自帶的SparseLU就夠用上萬自由度以上可以接PARDISO或MUMPS這類專業(yè)求解器。XFEM的矩陣天然有富集自由度帶來的條件數(shù)問題直接求解器比迭代法穩(wěn)得多省去調(diào)預(yù)條件子的時(shí)間。實(shí)測下來迭代法在富集節(jié)點(diǎn)分布不均時(shí)經(jīng)常跑不動(dòng)而直接求解器只是略慢但結(jié)果非常穩(wěn)。3.4 后處理與VTK可視化算完矩陣只能得到一堆自由度系數(shù)要把位移場、應(yīng)變場、裂縫形態(tài)可視化出來我一般直接輸出VTK格式交給ParaView或自己寫個(gè)小瀏覽器看。注意XFEM的后處理不能跟標(biāo)準(zhǔn)有限元一樣只把節(jié)點(diǎn)位移畫出來。因?yàn)楦患杂啥刃枰饐卧摹盎謴?fù)”到節(jié)點(diǎn)物理值上。最穩(wěn)妥的做法遍歷所有單元對每個(gè)單元的高斯點(diǎn)先判斷該點(diǎn)是否在裂縫影響區(qū)內(nèi)再根據(jù)水平集值是否跨越裂縫決定最終位移表達(dá)。輸出VTK時(shí)可以把富集單元內(nèi)部的位移剖成多個(gè)子像素點(diǎn)用等值線畫出來否則裂縫兩側(cè)的位移跳躍畫不出來視覺效果上就像裂縫不存在一樣。代碼里要定制一下數(shù)據(jù)導(dǎo)出把每個(gè)節(jié)點(diǎn)的水平集值一并輸出檢查水平集場是否跟預(yù)期一致。這一步價(jià)值巨大——有一次我發(fā)現(xiàn)某個(gè)節(jié)點(diǎn)的水平集符號(hào)反了就是靠把φ場導(dǎo)出來在可視化軟件里看出來的比對著數(shù)值查快得多。4. 復(fù)雜裂縫建模的實(shí)操技術(shù)4.1 多裂縫相交與分叉處理復(fù)雜裂縫網(wǎng)絡(luò)里裂縫相交、分叉是常態(tài)。很多初學(xué)XFEM的人會(huì)想一條裂縫我都處理明白了多條不就是復(fù)制粘貼嗎真做起來發(fā)現(xiàn)不是這么回事。兩條裂縫相交交點(diǎn)把裂縫一分為二交點(diǎn)上方的單元位移跳躍的定義變得微妙因?yàn)榱芽p一側(cè)的位移跳躍可能同時(shí)受兩條裂縫影響。工程上常用的折中方案是先做布爾運(yùn)算把相交點(diǎn)插入裂縫分段再對穿過同一單元的多個(gè)裂縫段做單元級(jí)“分割”把一個(gè)單元內(nèi)部按裂縫段切成若干子區(qū)域每個(gè)子區(qū)域單獨(dú)積分。這個(gè)“單元內(nèi)部再分片”的操作其實(shí)是XFEM框架下最重要的技術(shù)細(xì)節(jié)之一。我在實(shí)現(xiàn)時(shí)先算出裂縫段與單元邊的所有交點(diǎn)把單元多邊形重剖成多個(gè)凸多邊形子區(qū)域子區(qū)域內(nèi)積分點(diǎn)單獨(dú)生成再用子區(qū)域的形函數(shù)值加權(quán)匯總。這條思路實(shí)現(xiàn)起來有點(diǎn)繁瑣但寫一次能應(yīng)付任意復(fù)雜的裂縫相交場景。我踩過的坑是兩個(gè)裂縫段交點(diǎn)落在一個(gè)單元邊上而不是內(nèi)部時(shí)容易漏掉一個(gè)子區(qū)域?qū)е履莻€(gè)單元沒被積分結(jié)果出現(xiàn)“裂縫經(jīng)過但位移場沒變化”的詭異現(xiàn)象。解決方法是做嚴(yán)格的多邊形裁剪算法測試把裂縫穿過單元的六種基本拓?fù)洳糠诌叴┻^、頂點(diǎn)穿過、邊到邊等全部設(shè)成用例跑一遍。4.2 裂尖單元的數(shù)值積分方案數(shù)值積分是XFEM里最考驗(yàn)?zāi)托牡牟糠?。被Heaviside富集的單元如果直接用高斯積分裂縫面的跳躍會(huì)讓積分產(chǎn)生振蕩因?yàn)楸环e函數(shù)在裂縫處不連續(xù)。工程上的標(biāo)準(zhǔn)做法是把該單元按裂縫面切分成兩個(gè)子區(qū)域每個(gè)子區(qū)域各自做高斯積分。對裂尖單元除了切分還要考慮分支函數(shù)在裂尖的√r奇異性影響——通常把含裂尖的子區(qū)域進(jìn)一步細(xì)分。實(shí)際編碼時(shí)我用的積分子程序流程是先判斷單元的富集類型普通單元直接2×2或3×3高斯積分Heaviside富集單元先跟裂縫段求交切成子多邊形再對每個(gè)子多邊形做T型/三角形剖分裂尖富集單元沿裂尖尖端區(qū)域額外細(xì)化。注意“子多邊形剖分”只用在高斯積分點(diǎn)的生成上不改變單元的拓?fù)渥杂啥?。也就是說切分純粹是數(shù)值積分的輔助手段切完的子片不引入任何新節(jié)點(diǎn)或新自由度。理解這一點(diǎn)很重要它決定了整個(gè)程序架構(gòu)的簡潔性。4.3 網(wǎng)格準(zhǔn)備與幾何交互雖然XFEM允許裂縫不跟網(wǎng)格走網(wǎng)格準(zhǔn)備也不是完全沒講究。裂縫和網(wǎng)格的相對位置決定了富集區(qū)域范圍也影響條件數(shù)。我的經(jīng)驗(yàn)是裂縫盡量避開節(jié)點(diǎn)特別是裂尖不要正好落在節(jié)點(diǎn)上否則水平集函數(shù)出現(xiàn)零值模糊富集處理容易出幺蛾子。裂尖落在單元內(nèi)部是最舒服的狀態(tài)裂尖富集節(jié)點(diǎn)分布在周圍一圈也均勻。另外一個(gè)經(jīng)驗(yàn)是裂縫跟單元邊的交角不要太刁鉆比如接近0度。如果一條裂縫幾乎貼著單元邊卻從單元內(nèi)部穿過切分出來的子區(qū)域會(huì)極其狹長積分誤差劇增。解決辦法是生成網(wǎng)格時(shí)有意控制單元尺寸和方向或者對裂縫幾何做一點(diǎn)點(diǎn)擾動(dòng)讓交角避開極端情況。這種網(wǎng)格“微調(diào)”不要求裂縫貼著網(wǎng)格但要求裂縫避免笨蛋式的共線調(diào)完好處非常明顯。模型邊界與裂縫相交也有講究。如果裂縫一頭伸到模型邊界上Heaviside富集就只在裂縫內(nèi)部的端點(diǎn)處成立邊界側(cè)的富集節(jié)點(diǎn)不能帶到邊界條件里否則會(huì)引入額外的反力路徑。處理方式是對裂縫與邊界交點(diǎn)的兩側(cè)節(jié)點(diǎn)單獨(dú)標(biāo)記出來不參與邊界約束裝配。這個(gè)問題我在做壓裂縫穿過井筒附近的案例時(shí)踩過一次排查了兩天。5. 滲透率處理從位移場到流動(dòng)響應(yīng)5.1 提取裂縫開度位移跳躍應(yīng)力場求解完成后最自然的產(chǎn)物就是每個(gè)裂縫段的位移跳躍值。所謂開度w就是裂縫面兩側(cè)位移向量的法向分量之差w(x) [u(x)] · n (u? - u?) · n實(shí)現(xiàn)時(shí)取裂縫段上若干采樣點(diǎn)利用富集自由度的解計(jì)算采樣點(diǎn)兩側(cè)的位移差。這里有一個(gè)工程中經(jīng)常遺忘的細(xì)節(jié)XFEM算出的是彈性位移場但裂縫實(shí)際開度還要減去初始開度。很多天然裂縫在無載荷狀態(tài)下已經(jīng)有一定張開的寬度作為永久的初始條件記錄在裂縫結(jié)構(gòu)體里應(yīng)力和應(yīng)變引起的增量開度疊加到初始開度上才能進(jìn)滲透率計(jì)算。如果裂縫面粗糙局部接觸還會(huì)導(dǎo)致開度局部為零甚至出現(xiàn)負(fù)開度。幾何上不允許裂縫面彼此嵌入所以在更新滲透率之前要對開度做截?cái)嗵幚硇∮谀抽撝档拈_度按殘余開度處理或者引入接觸判定。這個(gè)截?cái)嚅撝翟趺炊ㄖ苯佑绊憹B透率結(jié)果苗頭很微妙。我做過一個(gè)參數(shù)敏感性分析同一模型閾值從1μm調(diào)到10μm等效滲透率能差出30%所以閾值要跟工程實(shí)際對照來選不要拍腦袋。5.2 立方定律與離散裂縫網(wǎng)絡(luò)滲透率張量裂縫開度到滲透率的橋梁是經(jīng)典的立方定律。它的假設(shè)是裂縫面為光滑平行板流體在其中做層流等效滲透率k_f w2 / 12單位是m2。這個(gè)公式先簡單又令人心生敬意。它成立的前提是裂縫面無限長、兩側(cè)完全平行、流動(dòng)完全層流。天然裂縫面有粗糙度壓裂裂縫面凹凸不平實(shí)際導(dǎo)流能力往往小于立方定律預(yù)測值。工程上常用折減系數(shù)把立方定律修正一下k_f (w2 / 12) · f_reduction折減系數(shù)f_reduction取值在0.1到1.0之間視粗糙度、填充物而定。要得到整個(gè)裂縫網(wǎng)絡(luò)的等效滲透率張量K對每條裂縫段做累積。二維情形下一條裂縫方向與整體坐標(biāo)系成θ角法向單位向量為n滲透率張量貢獻(xiàn)為K_f (w3 / 12) · (δ - n ? n) / L這里w的三次方很要命它意味著大開度的裂縫主導(dǎo)整個(gè)網(wǎng)絡(luò)的滲透率也意味著開度估算誤差會(huì)以三次方的倍率放大到滲透率里。所以前面講的初始開度、粗糙度修正、位移跳躍提取精度每一樣都馬虎不得。我調(diào)試時(shí)試過一些裂縫分布某條大裂縫開度差0.1mm等效滲透率就翻倍。想讓結(jié)果穩(wěn)定開度提取那塊需要做局部應(yīng)力場的平均或光滑化不要直接用單點(diǎn)采樣值。5.3 流固耦合的迭代思路搞完滲透率張量接著就要算流動(dòng)響應(yīng)。只算單次應(yīng)力載荷下的滲透率那是“靜態(tài)映射”真實(shí)項(xiàng)目往往需要流固耦合孔隙壓力變化改變有效應(yīng)力有效應(yīng)力改變裂縫開度開度再改滲透率如此往復(fù)。耦合方式有兩種。全耦合monolithic把固體位移和流體壓力放進(jìn)同一個(gè)非線性系統(tǒng)里矩陣維度大、收斂難度高。交錯(cuò)迭代staggered更友好先固定壓力場求解位移場再固定位移場求解壓力場循環(huán)直到收斂。對工程分析交錯(cuò)迭代已經(jīng)夠用而且程序上是在已有的XFEM固體求解器和滲流求解器之間加一個(gè)數(shù)據(jù)轉(zhuǎn)接層比推倒重來簡單太多。我用交錯(cuò)迭代時(shí)發(fā)現(xiàn)剛度比裂縫導(dǎo)流能力vs基質(zhì)滲透率特別懸殊時(shí)收斂會(huì)變慢甚至振蕩。加一個(gè)亞松弛因子比如新開度 0.7×新開度 0.3×舊開度能讓收斂穩(wěn)定不少。亞松弛因子選多大沒有標(biāo)準(zhǔn)答案跟問題有關(guān)我一般從0.5開始試觀察迭代殘差曲線逐漸加大到收斂最快且不振蕩的臨界值。6. 我踩過的坑常見問題與調(diào)試實(shí)錄6.1 富集自由度編號(hào)錯(cuò)亂導(dǎo)致矩陣奇異第一個(gè)大坑出現(xiàn)在自由度裝配上。表面上看邏輯很簡單但一旦裂縫穿過多個(gè)單元同一個(gè)節(jié)點(diǎn)可能同時(shí)被Heaviside富集和裂尖富集標(biāo)記編號(hào)循環(huán)里就要保證兩種富集的自由度都分配好。我第一次調(diào)試時(shí)裂尖富集節(jié)點(diǎn)在循環(huán)末尾被重復(fù)編號(hào)導(dǎo)致一個(gè)富集自由度指向了另一個(gè)節(jié)點(diǎn)的自由度整個(gè)剛度矩陣在未施加約束前就有了一列相同的列結(jié)果求解器報(bào)“矩陣奇異”。排查方法沒什么捷徑寫一個(gè)單元級(jí)dof映射表檢查腳本逐個(gè)單元對比它引用的自由度編號(hào)是否在全局自由度范圍內(nèi)且沒有多個(gè)單元使用不同含義的同一編號(hào)。后來我把編號(hào)邏輯提取成一個(gè)獨(dú)立函數(shù)任何富集類型變化時(shí)先跑一遍單元連線測試能提前抓到90%的錯(cuò)亂問題。6.2 數(shù)值積分不當(dāng)造成的應(yīng)力震蕩第二個(gè)坑是應(yīng)力震蕩。裂縫穿過的單元?jiǎng)傞_始我圖省事沒有做單元子區(qū)域剖分想當(dāng)然用高密度高斯點(diǎn)硬剛。結(jié)果裂縫面兩側(cè)的應(yīng)力場出來了自己都沒見過的鋸齒狀分布尤其在裂尖附近應(yīng)力等值線跟得了瘧疾似的瘋狂抖動(dòng)。這背后的道理其實(shí)簡單Heaviside函數(shù)在裂縫面上跳躍而被積函數(shù)在某個(gè)高斯點(diǎn)處有階躍高斯求積公式對不可導(dǎo)的間斷函數(shù)的逼近力不從心。要做單元級(jí)“內(nèi)部劃片切割”把裂縫面兩側(cè)的高斯積分點(diǎn)徹底分開讓積分點(diǎn)永遠(yuǎn)不跨過間斷線。做完這個(gè)改進(jìn)應(yīng)力分布肉眼可見地平滑了。我建議判據(jù)是看裂尖單元的K_I值在相鄰幾個(gè)單元間的連續(xù)性如果跳來跳去優(yōu)先懷疑積分方案。還有個(gè)細(xì)節(jié)裂尖分支函數(shù)在θ±π處有間斷atan2的返回值決定了θ的正負(fù)裂尖富集單元的積分子區(qū)域必須按這個(gè)間斷線再切開。這個(gè)坑非常隱蔽因?yàn)榱鸭馓幬灰茍霰旧磉B續(xù)你只觀察位移場看不到問題但應(yīng)變場會(huì)無比混亂。6.3 滲透率張量與邊界條件不匹配第三個(gè)坑來自滲透率張量的驗(yàn)證。用等效滲透率張量去算宏觀流動(dòng)結(jié)果跟解析解差距離譜。一開始我以為平方定律的折減系數(shù)設(shè)錯(cuò)了后來發(fā)現(xiàn)是張量的坐標(biāo)變換出了問題。裂縫的網(wǎng)絡(luò)方向定義在全局坐標(biāo)系但滲流求解器的邊界條件通常施加在模型邊界上兩者并不總是一致——問題不在物理模型而在坐標(biāo)系的變換處理。解決方法是寫一段獨(dú)立的驗(yàn)證用例單條裂縫橫穿方形模型已知其傾角、開度手算等效滲透率張量再讓程序跑一遍對比。這種“元測試”看起來簡單直接卻能一下子抓住變換矩陣的符號(hào)錯(cuò)誤、索引錯(cuò)位等低級(jí)問題。從那以后我每改一次數(shù)據(jù)結(jié)構(gòu)都會(huì)把這些最小用例重新跑一遍確認(rèn)沒有回歸。6.4 后處理可視化裂縫“消失”了的假象最后一個(gè)是后處理坑。輸出VTK文件后在桌面軟件里打開裂縫面明明穿過了好幾個(gè)單元但等值線圖上裂縫線消失了。當(dāng)時(shí)我還以為是富集項(xiàng)沒生效后來發(fā)現(xiàn)是后處理腳本里只讀取了標(biāo)準(zhǔn)節(jié)點(diǎn)位移根本沒把富集自由度加回去。這里提醒大家XFEM后處理絕不能套標(biāo)準(zhǔn)有限元的顯示流程。你要為每個(gè)單元生成一套“顯示用”的局部細(xì)化網(wǎng)格逐點(diǎn)把富集項(xiàng)拼進(jìn)位移值里。以二次元一個(gè)等值線剖視圖為例先把該單元內(nèi)部按虛擬細(xì)分成N×N個(gè)子格每個(gè)子格頂點(diǎn)都判斷水平集符號(hào)然后根據(jù)符號(hào)把對應(yīng)的富集位移加進(jìn)去這樣畫出來的圖里裂縫面兩側(cè)的位移跳躍才肉眼可見。最后分享點(diǎn)個(gè)人體會(huì)整個(gè)項(xiàng)目做完我最大的感受是XFEM不是“換一套網(wǎng)格”而是“換一套數(shù)據(jù)組織方式”。一旦你把幾何水平集、場位移富集、物理開度、滲透率分開管理復(fù)雜裂縫再難也只是數(shù)據(jù)結(jié)構(gòu)的復(fù)雜性而不是算法上的不可逾越。我自己在編碼過程中最大的收獲是養(yǎng)成了“任何富集操作都要配一個(gè)最小驗(yàn)證用例”的習(xí)慣——這個(gè)習(xí)慣看起來笨拙但它讓我避開了數(shù)不清的隱藏bug。如果你正準(zhǔn)備上手這類項(xiàng)目我建議先別急著追求復(fù)雜的多裂縫相交案例先做一條單裂縫、一個(gè)規(guī)則網(wǎng)格、一個(gè)已知解析解驗(yàn)證跑通之后再逐步加復(fù)雜度。相信我這一套走完再看復(fù)雜裂縫網(wǎng)絡(luò)會(huì)輕松很多。