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