COMSOL與AVIZO聯合仿真土體三維細觀滲流特性

蔡沛辰， 闕云， 李顯

(福州大學土木工程學院， 福建 福州 350108)

0 引言

隨著計算機技術的迅速發展， 不同領域多學科信息交叉的聯系越來越緊密， 僅靠單一領域已不足以解決復雜問題， 如流固耦合、 邊坡失穩變形等[1-2]. COMSOL Multiphysics 是一款多物理場耦合的數值仿真軟件， 具有靈活、 易用、 易于擴展的特性[3]. 目前， 與其聯合仿真的軟件主要包括MATLAB、 SOLIDWORKS、 CAD三種. 如： 1) 基于MATLAB軟件. 金亮等[4]提出將MATLAB與COMSOL聯合進行永磁發電機的優化仿真， 并實現兩者之間數據的傳遞； 李晶晶[5]通過COMSOL與MATLAB之間相互調用， 構建水下腔體的磁化特性仿真系統， 并對結果進行誤差分析. 2) 基于SOLIDWORKS軟件. 范怡哲[6]聯合COMSOL有限元軟件對軌道交通受電弓三維模型進行力學特性及故障分析. 3) 基于CAD軟件. 褚靜等[7]聯合COMSOL軟件對LED燈珠進行熱能分布的研究. 上述相關研究已較為成熟， 但目前SOLIDWORKS、 CAD都只局限于構建人為設定的幾何模型， 對于自然界中不可準確定量表征的物體， 實現存在較大困難， 如真實土體孔隙結構模型. 此外， 應用較為廣泛的MATLAB軟件， 雖可采用代碼編程來構建復雜模型， 但無法直觀地對計算模型進行可視化處理和網格修復等操作. 因此， 本文提出COMSOL與AVIZO聯合仿真方案， 對真實土體的三維模型仿真過程進行詳細的可視化研究. AVIZO是針對地質和材料科學的一款強大可視化軟件， 具有高級三維可視化、 圖像處理、 三維重構、 網格劃分及修復等功能[8]. Bird等[9]、 Fang等[10]、 王平全等[11]和苗杰[12]采用此方法研究過巖石、 煤體的滲透率問題， 但對于聯合仿真具體實現過程未詳細描述， 不利于其他學者參考學習， 同時對土體領域的三維模型聯合仿真研究鮮有提及.

鑒于此， 本研究以原狀花崗巖殘積土為對象， 探討COMSOL與AVIZO聯合仿真技術的具體實現過程， 并以構建的三維土體模型滲流場仿真為例， 驗證聯合仿真方法的可行性， 為現有土體三維重構仿真模擬及流體滲流研究提供新途徑.

1 CT掃描及三維模型構建

1.1 CT掃描及圖像處理

試驗原狀土選自福州市某地山坡， 現場取樣如圖1所示. 選取植被茂密的位置， 取樣前先用鐵鍬清理表面腐殖層， 清理面積為100 cm×100 cm， 厚度為20 cm， 最終獲取尺寸為15 cm×15 cm×40 cm的土柱試樣， 試樣的基本物理參數如表1所示. 對獲取的土柱試樣進行工業CT掃描試驗， CT掃描試驗設備名稱為C450KV高能量工業CT， 工作電壓為450 kV， 電流為63 mA， 掃描最低分辨率為0.15 mm， 掃描后得到一系列CT掃描圖像， 如圖2所示.

圖1 現場取樣Fig.1 Field sampling

圖2 CT掃描切片Fig.2 CT scanning slice

表1 所選試樣的基本物理參數

現實中， 圖像在成像或傳輸過程中， 常受到設備與外部環境等因素影響[13]， 故在構建合適的三維模型前需對圖像進行必要處理， 主要步驟包括： 1) 圖像去噪. 采用AVIZO中Unsharp masking和Sobel功能提高圖像對比度及檢測圖像邊緣， 同時對噪聲具有平滑和抑制作用. 2) 二值化閾值分割. AVIZO中Watershed Segmentation可根據圖像的灰度值不同， 將圖像分割成土壤基質和孔隙兩種區域， 圖3為閾值分割圖， 其中， 藍色為分割后的孔隙區域. 對比可知分割過程非常準確地捕捉到了孔隙結構的邊界， 為后續孔隙網絡模型的構建提供一定基礎.

圖3 閾值分割Fig.3 Threshold segmentation

1.2 三維模型構建

在上述圖像處理基礎上， 以孔隙居中軸線體系為基礎[14]， 通過膨脹算法進行孔隙的分割和提取[15]， 最終構建成孔隙網絡模型. 圖4為AVIZO中孔隙網絡模型構建界面圖， 其中左側為操作流程和參數設置界面， 中間為可視化窗口， 右側為孔隙參數表. 模型中球形代表孔隙結構， 連接球形的棍代表兩個孔隙之間的喉道， 對不同的孔隙大小采用不同的顏色進行渲染， 從圖4中可清楚地看出土體孔隙連通情況.

2 COMSOL與AVIZO可視化數據對接過程

采用AVIZO軟件中模型簡化—網格劃分—網格優化及修復—網格測試等操作來獲得高質量的三維模型網格， 并輸出可導入COMSOL的STL模型文件， 為后續聯合仿真提供良好的基礎. 此外， COMSOL對真實三維孔隙結構進行仿真時， 計算量大且需反復對網格進行修復， 操作過程極為復雜， 故截取模型中尺寸6 mm×6 mm×6 mm (40體素×40體素×40體素)的孔隙連通性較好區域進行模擬.

2.1 模型簡化

三維模型中各孔隙結構極為復雜， 若直接將其導入COMSOL中進行網格劃分， 將會導致生成的網格質量較差或者網格劃分出現錯誤， 進而造成后續的仿真模擬無法實現， 因此在仿真前需對模型進行簡化處理.

圖5為三維孔隙模型中的孤立點狀孔隙分布圖. 從圖5中可發現： 三維模型中存在部分孤立的點狀孔隙， 這些孔隙孔徑小、 數量多， 且對土壤材料的滲透性能影響較小， 加之STL文件可能存在數據冗余現象. 故為減少后續處理過程中的計算量， 在網格劃分前需要通過孤立孔隙的刪除功能來清除這些孤立的點狀孔隙. 具體操作過程為： Segmentation→Remove islands→3D volume→Highlight all islands→Apply. 點狀孤立孔隙刪除后， 如圖6所示.

圖5 孤立點狀孔隙分布圖Fig.5 Isolated dotted pore distribution map

圖6 孤立點狀孔隙刪除示意圖 Fig.6 Schematic diagram of isolated dotted pore deletion

2.2 網格劃分

經前處理簡化后， 采用編輯器即可完成三維模型的網格劃分， AVIZO軟件進行網格劃分時， 默認在平坦區域生成的網格稀疏， 而對于曲率較高的區域， 生成的網格較密集， 其作用是： 既可以減少網格數量， 又能較好地保留三維模型的原始結構. 圖7為土壤基質和孔隙結構的網格劃分. 網格劃分操作過程為： Simplification Editor→Generate Tetra Grid→Max dist→Simplify now→Close. 需要注意的是： 1) Simplify now功能可能會丟失部分精確的初始表面， 因此在簡化之前需復制原表面. 2) Max dist最優取值為模型實際最小體積尺寸的1%左右， 體積的實際尺寸可以使用Local Axes來查看.

圖7 網格劃分圖Fig.7 Meshing diagram

2.3 網格優化與再修復

真實原狀土孔隙模型在網格劃分過程中會生成一些質量較差的網格， AVIZO提供了網格檢測與修復的功能. 將劃分好的網格運行寬高比測試， 測試模型網格中寬高比不合理的三角形， 并顯示如圖8所示紅色標識. 檢測過程具體為： Surface Editor→Surface menu→Tests submenu→Intersection test.

圖8 測試出的劃分不合理網格圖Fig.8 Unreasonable grid graph

此外， AVIZO可對檢測出的問題網格進行優化修復， 提供平移頂點、 合并頂點等功能. 修復過程具體為： 1) 平移頂點. Surface Editor toolbar→Translate Vertices→Translate， 如圖9(a)所示. 2) 合并頂點. Surface Editor toolbar→Contract Edges→Remove， 如圖9(b)所示.

2.4 模型網格質量測試及數據對接

網格修復完成后， 通過AVIZO網格質量測試的三維模型可以保存為STL格式文件， 進而與COMSOL進行數據交互對接， 至此COMSOL與AVIZO三維模型可視化數據對接過程完成， 反之， 需重復上述步驟， 重新對網格進行劃分.

3 COMSOL與AVIZO可視化聯合仿真應用實例

為驗證COMSOL與AVIZO聯合仿真技術的可行性， 按照上述數據對接方法， 將導出的STL模型文件導入COMSOL軟件， 進行三維模型滲流模擬的聯合仿真， 其中假設孔隙中的流體為不可壓縮流體， 滲流狀態為層流. 圖10為經過數據對接導入COMSOL中的幾何模型及網格質量分布圖， 從圖10中可見： 模型網格劃分質量分布較好， 大都集中分布于0.7以上.

圖10 幾何模型及網格質量分布Fig.10 Geometric model and mesh quality distribution

3.1 控制方程

質量和動量傳遞可采用N-S斯托克斯方程作為理論基礎[16]， 控制方程如下:

(1)

式中：ρw為流體密度， kg·m-3；p為壓力， Pa；I為單位矩陣；μ為流體動力粘度， Pa·s；F為體積力， N·m-3；u為流速， m·s-1；g為重力加速度， m·s-2.

3.2 基礎參數設定

采用水作為滲流流體， 并在滲流模塊中設置基本參數， 具體參數如表2所示.

表2 邊界條件參數

3.3 邊界條件

以Z方向滲流模擬為例， 進行邊界條件設定： 沿Z軸， 兩個相對面分別設置為入口和出口， 4個側面

為自由滑移壁， 其余內壁為無滑移狀態， 如圖11所示.

圖11 邊界條件Fig.11 Boundary condition

在入口和出口處采用水壓邊界條件， 即：

P=Pin，n·μ▽2u=0

(2)

P=Pout，n·μ▽2u=0

(3)

水壓加載邊界處：

u=0

(4)

側面邊界：

u·n=0

(5)

t·μ(-▽P+▽2u)n=0

(6)

式中：Pin、Pout分別為入口和出口水壓值， Pa；n為方向向量；t為滲流時間， s.

3.4 模擬結果及分析

選擇瞬態求解器， 求解時間(0, 0.1, 2 s)， 計算自由度個數 526 096(加1個內部自由度)， 采用后處理模塊中三維繪圖組繪制滲流過程中流線分布圖和壓力場分布圖， 如圖12和圖13所示. 從圖12中可看出： 速度分布流線并沒有覆蓋整個模型空間， 分析原因是孔隙幾何形狀的復雜性會不可避免地導致某些孔隙凝滯， 即在這些孔隙中流體未滲入. 從圖13可發現： 沿滲流方向， 孔隙所受的壓力逐漸降低， 且孔徑越小， 所承受壓力越大.

圖12 不同方向流線分布圖Fig.12 Streamline distribution diagram in different directions

圖13 不同方向壓力場分布圖Fig.13 Distribution map of pressure field in different directions

滲透率定義為一定壓差下， 多孔介質允許流體通過的能力， 計算公式如下式所列[12]：

(7)

式中：L為土樣流體流動方向長度， m；A為土樣橫截面積， m2； Δp為土樣壓降， Pa.

最后， 以Z方向滲流模擬為例， 對出口邊界滲流速度進行體積分， 可得到通過土樣的體積流量Q為6.08×10-12m3·s-1， 再由式(7)可求得Z方向的滲透率K的結果為1.03×10-6μm2， 相比文獻[17]中對原狀花崗巖殘積土滲透率測試結果2.43×10-6μm2(壓差1 kPa)， 聯合仿真技術測得滲透率K略偏小， 但仍處于同一數量級， 同時參考文獻[18]， 可進一步驗證COMSOL與AVIZO聯合仿真結果的正確性.

4 結語

基于CT掃描與圖像處理技術， 構建真實土體的3D模型， 并依據COMSOL與AVIZO數據交互對接技術對原狀土單相水滲流過程進行模擬. 研究發現， 可視化聯合仿真技術能較好地模擬原狀土體三維模型的單相水滲流， 同時驗證了兩者數據交互對接結果的正確性. 該研究可為今后土體三維重構、 細觀孔隙結構定量化表征及流體滲流研究提供一定借鑒.