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