數值模擬顯式計算方法

『壹』 數值模擬模型

6.2.2.1 地下水系統的模擬模型

潛水含水層的模擬模型：

含有協變數的地下水動態規劃管理模型研究

式中：D——地下水系統的模擬滲流區域；

h——地下水位（m）；

M——含水層厚度（m）；

K——滲透系數（m/d）；

T——導水系數（m²/d）；

t——時間（d）；

μ——給水度；

μ^*——承壓含水層貯水系數；

Γ₁，Γ₂——一類及二類邊界；

n→——邊界的外法線方向；

h₁（x，y，t）——第一類邊界的地下水位（m）；

h₀（x，y）——初始地下水位（m）；

q（x，y，t）——第二類邊界的單寬流量（m²/d）；

P——地下水開采量（m³/d）；

ε——天然條件的源匯項（m³/d），如降水入滲補給、蒸發排泄等；

Q——與地下水位有關的源匯項流量（m³/d），即協變數。

6.2.2.2 模擬軟體選取

本次工作採用Visual MODFLOW軟體對水流進行模擬，該軟體具有處理含有協變數的地下水模擬模型的功能。

6.2.2.2.1 河水與地下水的交換量

Visual MODFLOW河流程序包輸入文件要求河流邊界條件的每個網格單元具有以下資料：

（1）河流水位：地表水體的自由水面標高，這個標高可以隨時間變化。

（2）河床底板：地表水體的滲流層（夾層）的層底板標高。

（3）水力傳導系數：一個表示由滲流層（河床）引起的地表水體和地下水間水流阻力的數值參數。

水力傳導系數（C）可以通過單元的長度（L）、單元中河流的寬度（W）、河床厚度（M）以及河床物質的滲透系數（K），利用下面公式計算得出：

含有協變數的地下水動態規劃管理模型研究

6.2.2.2.2 排水溝的排泄量

Visual MODFLOW支持MODFLOW包括的標准排水溝邊界程序包。排水量正比於含水層的水頭與某一固定水頭或高程之差。如果含水層水頭低於排水溝的固定水頭，排水溝程序包假設排水就沒有效果。

排水溝程序包要求包含這個邊界條件的每個單元具備下列輸入資料：

（1）水頭標高：排水溝標高或者排水溝內的自由水面的排水溝水頭。假設排水溝並沒有充滿水，這樣排水溝管內的水頭接近等於排水溝管的中間標高。

（2）水力傳導系數：排水溝水力傳導系數是一個綜合系數，描述了排水溝和地下水系統之間的水頭損失。計算排水溝的水力傳導系數沒有通用的公式，大多數情況下，計算排水溝水力傳導系數所需的詳細資料對於地下水模型來說無法得到。這些細節包括排水溝周圍水頭分布的詳盡資料、排水溝附近含水層的滲透系數、填充材料的分配、排水管的數目和尺寸大小、堵塞物質的總量和堵塞物質的滲透系數。通常利用水流率測量值和水頭差來計算排水溝傳導系數。排水溝水力傳導系數在模型校正時通常會做調整。

6.2.2.2.3 蒸發量

蒸發蒸騰程序包模擬植物蒸騰、直接蒸發以及由地下水飽水狀態移走水而引起的地表面滲漏作用。

蒸發蒸騰程序包需要以下信息：

（1）蒸發蒸騰速度：當水面高程與網格單元高程一致時的蒸發蒸騰速度。

（2）極限深度：網格單元頂板以下的深度，在此深度蒸發蒸騰速度可以忽略不計。

蒸發蒸騰程序包方法基於以下假設：

（1）當水面位於或者高於地表面（第一層頂板）時，水面蒸發蒸騰損失以用戶定義的最大量進行。

（2）當水面高程低於極限深度時，或者位於第一層以下時，水面蒸發蒸騰可以忽略不計。

（3）兩個極限之間，水面蒸發蒸騰隨水位高程呈線性變化。

6.2.2.2.4 泉流量

Visual MODFLOW中沒有單獨計算泉流量程序包，泉流量在識別、檢驗和預報階段均採用排水溝計算模塊計算。

應用排水溝程序包計算泉流量時要求輸入下列資料：

（1）水頭標高：輸入泉出露點地表高程。

（2）水力傳導系數：輸入經識別、檢驗後與排水溝水力傳導系數等效的泉流量的比例系數。

在本次研究過程中，對松原地區地下水系統中的互饋協變關系問題作了一些合理的概化，重點考慮了泉流量、蒸發量、地下水與河水的交換量這三個協變數，並採用Visual MODFLOW軟體相應的程序包計算。

6.2.2.3 研究區網格剖分

對松原地區水文地質圖掃描，並應用ArcGIS軟體矢量化，作為計算區的剖分底圖。計算模擬區應用Visual MODFLOW軟體共剖分出80行、100列，4071×3個矩形單元體（圖6.6）。

圖6.6 研究區網格剖分圖

6.2.2.4 模型的校正

建立一個描述實際地下水系統的數值模擬模型，實際上就是找一個描述它的偏微分方程並確定其定解條件。一個正確可靠的數學模擬模型應當是實際地下水系統的復製品。也就是說，當施加天然或人為影響時，模擬模型的反應與實際地下水系統的反應應當非常接近，但在實際工作中很難直接達到這一點。首先，選用的方程的類型不一定合適。其次，代入方程中的各種參數，不論從現場還是在實驗室都難以准確的獲得，邊界條件也往往缺乏可靠的資料。此外，在建立水文地質概念模型的過程中，對地下水系統的實際條件又做了許多的簡化和假設。所以上述這些因素都可以導致初步建立的模擬模型與實際問題有很大的差距。因此，我們必須根據抽水試驗和已有的地下水系統長期動態觀測資料對初步建立的模型進行校正。

模型的校正包括兩個過程，即識別和檢驗，通常運用已有的不同時段的兩套地下水數據，一套用於識別，另一套用於檢驗。然後，把數值模擬模型計算出的結果與已有的實際數據進行對比，看兩者是否一致。若不一致，就要對模型繼續校正，直到滿意為止。

6.2.2.4.1 模型的識別

（1）模型識別時段的選擇：本次模擬的識別時段選取1999年12月1日到2000年2月29日，共計91d。本次模擬選取枯水期，因為期間源匯項少，地下水位比較穩定，較容易反映含水層的水文地質特性。

（2）參數分區與初值的確定：根據區域水文地質情況和抽水試驗成果，潛水含水層水文地質參數分為9個區，弱透水層水文地質參數分為4個區，承壓含水層水文地質參數分為8個區。參數分區見圖6.7～圖6.9。

圖6.7 潛水含水層參數分區圖

圖6.8 弱透水層參數分區圖

圖6.9 承壓含水層參數分區圖

（3）源匯項的計算：模型的識別時段正處於冬季，所以模擬區內無大氣降水入滲補給、灌溉回滲補給、農業開采和蒸發排泄等，排泄方式主要有開采排泄和泉的排泄，源匯項的計算較簡單。

（4）模型識別時段擬合結果：從水位擬合圖（圖6.10～圖6.11）可以看出，潛水觀測井水位誤差的平均絕對值為0.295m，擬合誤差小於0.5m的觀測井占總觀測井的88.73%；承壓水觀測井水位誤差的平均絕對值為0.302m，擬合誤差小於0.5m的觀測井占總觀測井的86.48%，可見模型識別取得了較好的結果。

表6.3—表6.5為潛水含水層、弱透水層和承壓含水層水文地質參數識別結果。

圖6.10 潛水含水層識別時段水位擬合圖

圖6.11 承壓含水層識別時段水位擬合圖

表6.3 潛水含水層水文地質參數識別結果表

表6.4 弱透水層水文地質參數識別結果表

表6.5 承壓含水層水文地質參數識別結果表

6.2.2.4.2 模型的檢驗

為了進一步驗證所建立的數值模擬模型和確定的水文地質參數的真實性，再次利用已有的地下水位動態觀測資料對模型進行檢驗。

（1）檢驗時段的選擇：選擇2005年6月1日至2005年8月31日，共計92d作為模型的檢驗階段。潛水含水層和承壓水含水層分別選取有長期觀測資料的觀測井進行模型的擬合檢驗。

（2）源匯項的計算：模型的檢驗時段正處夏季豐水期，源匯項比較復雜。主要補給來源有大氣降水入滲補給、灌溉入滲補給以及側向徑流補給等。主要排泄方式有人工開采排泄、潛水蒸發排泄以及泉的排泄。

（3）模型的檢驗時段擬合結果：從模型的檢驗時段擬合結果可以看出（圖6.12，圖6.13），潛水含水層水位擬合誤差小於0.5m的觀測井占總觀測井的89.53%；承壓含水層水位擬合誤差小於0.5m的觀測井占總觀測井的85.28%。模型的檢驗結果表明，建立的地下水數值模擬模型和確定的水文地質參數是真實可靠的。

6.2.2.5 模型的預報

經過識別和檢驗的模型，就相當於實際地下水系統的復製品，就可以以模型的行為代表實際地下水系統的行為了。對於各種自然因素、人為因素的影響（比如各種規劃開采方案），都可以通過數學模擬模型的運轉來預報其後果，這也是建立模擬模型的主要目的之所在。

6.2.2.5.1 預測期、初始條件的確定以及邊界條件的處理

本次預測期限為4年，即2008～2012年。根據已有的地下水位動態觀測資料，選擇地下水位預報的初始時刻為2008年9月1日。

圖6.12 潛水含水層檢驗時段水位擬合圖

圖6.13 承壓含水層檢驗時段水位擬合圖

由於開采水源地距離計算模擬區邊界較遠，並且考慮到開采時對邊界的影響較小，因此本次採用初始時刻的邊界條件進行計算。

6.2.2.5.2 源匯項的預報

（1）降水量和水面蒸發量的預報：本次應用蒙特卡羅（Monte-Carlo）方法對研究區內的降水量和水面蒸發量進行隨機模擬和預報，這是一種應用隨機數來進行模擬試驗的方法。該方法把降水量和蒸發量作為一個隨機變數來處理，根據過去的資料提取研究區降水量和蒸發量所具有的統計規律性，這些統計規律性代表了該地區降水量和蒸發量的本質特徵。

（2）灌溉入滲補給量的預報：灌溉入滲補給量根據研究區內的農業用水額度來進行預報。

（3）開采量的預報：由於經濟發展和人口增長等原因，松原地區用水量在逐年增加。2002年地下水開采量為4.×10⁸m³ ，2008年地下水開采量為5.2201×10⁸m³ ，平均年增長率定為3.83%。按該增長率對2012年開采量進行預測，即6.067×10⁸m³。6.2.2.5.3 預報結果及分析

根據地下水開采量平均年增長率3.83%對開采量進行預測，即計算模擬區2012年地下水開采量為6.067×10⁸m³。從預報結果的水位等值線圖（圖6.14和圖6.15）可以看出，潛水水位和承壓水位都有一定程度的下降，部分地區水位下降明顯。因此，對水位持續下降地段，應採取科學管理和一定的控制措施。此外，由於地下水開采量的增加以及當地混合開采地下水，致使潛水層和承壓水層直接產生了水力聯系，承壓含水層的承壓性質正在減弱，向潛水轉化。

圖6.14 2012年9月1日潛水等水位線圖

針對上述情況，在地下水開發利用中應採取一定的對策：①要科學地布井，有計劃地進行淺、中、深層分層開采地下水；②開展地下水長期監測，以便掌握區域水情的變化趨勢，發現問題及時解決；③加強水資源管理工作，合理開發利用地下水資源，確定開采量的大小和開採的時段長度，控制超量開采局面。

圖6.15 2012年9月1日承壓水等水位線圖

『貳』 數值模擬方法

如前所述，基於鏡質組化學反應動力學的煤熱演化史數值模擬方法經歷了簡單函數關系模式、受熱時間－經驗法模式、反應活化能－溫度函數模式和平行反應化學動力學模式四個發展階段（秦勇等，1995）。其中，目前在化石能源地質界最為常用的是TTI法、LOM法和Tissot法，近年建立起來的EASY%R_o方法也已成為數值模擬的重要發展方向。

（一）時間－溫度指數（TTI）法

這種方法的基本思想由Lopatin（1971）首先提出，後經Waples（1980）根據31個盆地402個樣品的實測資料加以完善，建立起了其數值模擬的方法體系，故又稱為Lopatin-Waples法。TTI法遵循化學反應動力學的基本法則來衡量溫度和時間兩個因素對煤化作用的貢獻，即煤化作用速度隨受熱溫度的增高而增大，溫度每增加10℃煤化反應速度增加一倍。

由此，可將TTI數學模式用連續函數加以表示：

山西南部煤化作用及其古地熱系統：兼論煤化作用的控氣地質機理

T=T_s＋G·D

式中：T——煤層受熱的古地溫溫度，℃；

T_s—－古地表溫度，℃；

G——古地溫梯度，℃/m;

D—煤層埋藏深度，m。

上述連續函數模式不便於常規運算，故在實際工作中通常採用分段積分方式，根據如下方程進行計算：

TTI＝∑γⁿ·△t

式中：γ——溫度效應因子，考慮到溫度每增高10℃煤級量值增加一倍的關系，取r=2;

n——溫度指數沖仿，其大小取決於溫度（表5-1）；

△t——煤層在某一溫度段的受熱時間，Ma。

表5-1TTI數學模式中溫度指數（n）與煤層受熱溫度段的對應關系

註：溫度每增（減）10℃，溫度指數增（減）1。

TTI值本身可以作為煤級或有機質成熟度指標，Waples（1980）也提供了TTI值與鏡質組最大油浸反射率的關系（表5-2）。此外，Wood（1988）通過回歸擬合，將二者之間關系採用下式表示：

lg%R_o＝—0.006（lgTTI）³+0.042（lgTTI）²+0.162（lgTTI）—0.397

表5-2時間溫度指數（TTI）與鏡質組反射率（R_o）的相關關系

應予指出，TTI法盡管幾十年來被化石能源地質界所廣泛採用，但仍存在某些明顯的不足有待於改進，國內外學者為此做過一定努力。例如，溫度每碧猜增加10℃煤化速率加快一倍的假設有其局限性，只能在活化能10～25kcal/mol（相當於20～160℃的溫度范圍，1kcal=4186.8J，下同）區間是可靠的（Magoon,1983）。再如，沉積有機質熱降解的反應活化能是受熱溫度的函數，隨成熟度的增加而逐步加大，不同類型有機質達到相同成熟度所需的活化能也有差異，而TTI模式將整個煤化過程中的反應活化能作為常量，顯然與實際情況有出入。

（二）悔判型有機成熟度水平（LOM）法

該方法由Hood（1975）建立，後經Bostick等（1979）完善而成為一種國內外廣泛採用的數值模擬方法，故又被稱為Hood-Bostick方法。

Hood等考慮到總體反應活化能隨受熱溫度而增高的變化趨勢（18～33kcal/mo1），採用鏡質組油浸反射率標定有機質成熟度水平（LOM）（圖5-1）。在模式中採用了有效受熱時間的概念，即溫度不低於最高受熱溫度15℃范圍內的受熱時間，建立起有機質成熟度、溫度及有效受熱時間之間的相互關系（圖5-2）。

圖5-1有機質成熟度水平（LOM）與其它有機質成熟度指標的關系（引自周中毅，1990）

Ⅰ—褐煤；Ⅱ₁—亞煙煤；Ⅱ₂—高揮發性煙煤；Ⅱ₃—中揮發性煙煤；Ⅱ₄—低揮發性煙煤；Ⅲ₁—半無煙煤；Ⅲ₂—無煙煤

1Btu=1055.06J

在實際工作中，LOM法採用圖解方法求取有機質成熟度或鏡質組反射率（圖5-2）。國內外應用效果表明，利用這種方法推算出的盆地古地溫溫度是比較可靠的，與由綜合研究所得出的結果以及盆地實際情況較為吻合（Vote,1981；周中毅等，1983,1984,1985）。

『叄』 數值模擬技術簡介

（一）研究現狀

地下多相、多組分流體運移數值模擬是在質量和能量守恆的基礎上，建立的多相流體運動以及反映地球化學運移擴散的數學模型，通過離散建立大量的線形或非線形方程組，然後利用計算機計算求解，再通過圖像顯示模擬結果，達到對工程問題和物理問題乃至相關其他問題研究的目的。CO₂地質封存數值模擬就是利用計算機模擬的方法，來解決CO₂進入地質封存系統後運移、轉化、水-岩-氣之間的相互反應、CO₂泄漏對淺部含水層影響及誘發的儲蓋層物性變化等一系列問題，從而指導CO₂地質封存工程的實施。

目前，國內外已開展的關於CO₂地質封存數值模擬的研究工作包括以下幾個方面：

1.超臨界CO₂-水多相流體運動模擬

Pruess等（2003）模擬了均質各向同性鹹水含水層中以恆定流量灌注CO₂條件下，灌注井井周非等溫徑向流情況。當忽略重力和慣性力效應時，模擬結果中存在相似變數ζ＝R²/t（其中，R為徑向流動距離，t為時間），CO₂飽和度、溶解CO₂質量分數、沉澱鹽的體積分數及流體壓力都是相似變數的函數。這與O' Sullivan（1981）及Doughty和Pruess（1992）的結果一致。兩相流的模擬考慮了CO₂和水的相對滲透率及毛細管力作用問題（Van Genuchten，1980），考慮了流體密度、黏度和CO₂溶解性隨壓力、溫度和含鹽度的變化，以及鹽的沉澱導致含水層滲透率的減小等因素。

Doughty和Pruess（2004）利用Fro鹹水含水層封存CO₂監測資料，反推了CO₂灌注後發生的物理和化學過程。他們採用TOUGH2數值模擬軟體對兩相（液、氣）三組分（CO₂、水和溶解NaCl）系統進行模擬。考慮15MPa和65℃條件下，超臨界CO₂在鹹水中為非混溶流體，並能部分溶解於鹹水的情況，分析了多相流系統邊界設定的影響及相對滲透率取值問題，即模擬中對側向邊界的設置為均開（或均閉），結果導致壓力的模擬結果與實際相比過低（或過高）。研究表明，由於上覆斷層對CO₂的封堵作用，側向邊界對CO₂彌散羽的影響不大。模擬結果還顯示，相對滲透率函數對CO₂彌散羽的演化有很強的影響。如何確定一個合適的相對滲透率以表徵CO₂灌注鹹水含水層的變化，仍是亟待解決的問題。Doughty和Pruess模擬了兩種CO₂殘余飽和度條件下的CO₂羽擴展隨時間的變化，發現存在較大差異。殘余飽和度較大的情況下，CO₂羽表現出緊縮狀，在浮力作用下運移較慢；相反，在殘余飽和度較小的情況下，CO₂羽流彌散很快，溶解性顯著提高。

2.多組分反應地球化學運移模擬

水-砂岩-CO₂相互作用往往形成一系列次生礦物，或次生礦物組合。Worden et al.（2006）通過岩石學和CO₂灌注長石砂岩的地球化學模擬表明，北海Magnus油田上侏羅統濁積亞長石砂岩中的鐵白雲石、高嶺石和石英可能具有成因聯系。其中，鐵白雲石中的碳來自有機成因的CO₂。Watson et al.（2004）通過CO₂氣與CH₄氣儲集砂岩的比較岩石學研究，證實澳大利亞Otway盆地Ladbroke Grove CO₂氣儲集砂岩中與CO₂氣灌注有關的次生礦物組合為鐵白雲石-高嶺石-次生石英。

Xu et al.（2005）利用一維砂岩-頁岩系統模型模擬了儲層中灌注的CO₂與礦物發生的化學反應過程，以及對儲層環境的影響。模擬顯示，在砂岩環境下，CO₂主要被方解石所固定，而方解石的沉澱導致孔隙度減小，進而導致滲透率相應減小。10萬年間，砂岩封存能力達到90kg/m³的封存能力，這些被礦物固定下來的CO₂可以永久封存。Zwingmann等運用地球化學模擬軟體EQ3/6進行的水-礦物-CO₂相互作用模擬也表明，若將CO₂灌注到日本本州島中北部新潟盆地更新世灰爪組砂岩，CO₂以溶於水和形成碳酸鹽礦物兩種形式封存，其中後者封存量最大為21.3mol/kgH₂O，可達總封存量的90％，形成的碳酸鹽礦物中也出現了片鈉鋁石。

3.耦合岩石力學模擬

從目前發表的論文及各國研究計劃的綜合報告上看，在CO₂鹹水含水層封存研究方面，對於CO₂運移機制的分析和模擬很少考慮應力場的耦合作用。事實上，CO₂灌注壓力和超臨界CO₂的浮力作用將改變地層應力狀態，即CO₂在上浮運移和側向擴散過程中，孔隙壓力可能會對原始裂隙和斷裂產生影響；CO₂在鹹水含水層中的長時期（千年級尺度以上）的封存，將改變含水層的地球化學狀態，CO₂-鹹水-含水層礦物的化學作用將可能導致岩體力學和水力學性質發生變化。

日本因位於4大板塊交界處與環太平洋構造帶中，活斷層密集發育，地震頻繁，地應力分布復雜，在CO₂地質封存評價方面，非常重視CO₂地質封存的力學穩定性研究（李琦等，2002；李小春等，2003）。李琦等（2002；2004；2006）提出了一個考慮初始地應力場、灌注壓力、CO₂浮力及含熱傳導作用的熱-水-力（THM）耦合模擬框架，考慮蓋層底部附近存在不同傾角斷層的二維平面應變地質封存問題。採用有限元演算法，對灌注CO₂流體對斷層穩定性的影響進行模擬分析。計算結果表明，為了避免斷層位移需要特別注意對灌注壓力的控制，因為CO₂灌注壓力對斷層滑動的影響遠大於CO₂羽流浮力帶來的影響。停止灌注CO₂後，CO₂羽流的上升則成為應力場擾動的主要因素。

（二）主要軟體介紹

近年來，計算機模擬技術在許多研究領域得到了廣泛的應用，開發出了許多優秀的模擬軟體和程序。同樣，可用於研究CO₂地質封存的數值模擬軟體也很多，主要有PHREEQC、GEM、ECLIPSE、TOUGHREACT、PetroMod、MUFTE-UG和NUFT等，它們都有各自的特點和適用性。在進行數值模擬之前，需對這些數值模擬軟體進行評價分析，選擇適用於所研究問題的模擬軟體。現對國際上常用的幾款軟體簡介如下。

1.PHREEQC

PHREEQC是一款用於計算多種低溫水文地球化學反應的計算機軟體。以離子締合水模型為基礎，PHREEQC可完成以下任務：(1)計算物質形成種類與礦物的溶解飽和指數；(2)模擬地球化學反演過程；(3)計算批反應與一維運移反應。另外，與多組分溶質-運移模型耦合的PHREEQC可生成PHAST，一個用於模擬地下水流系統的三維反應-運移模擬器。但由於PHREEQC是在單相水流的基礎上建立的模型，因而不能模擬超臨界CO₂-水的兩相流運動。

PHREEQC最簡單的應用就是計算溶液中各種化學物質的分布，以及溶液中礦物與氣體的飽和狀態。反演模擬功能可推導和量化在流動過程中，能夠反應化學物質變化的化學反應方程。PHREEQC可處理的反應方程包括建立礦物、表面配合物、陽離子交換劑、土壤溶液、氣體組分單位分壓、給定壓力或給定體積氣相間平衡的物質運移反應。在模擬這些均衡反應的同時，PHREEQC還可以模擬動力化學與生物反應，以及模擬從簡單的線性衰變（代謝物降解或放射性衰變）到復雜的依賴於溶液化學組成和微生物數量確定的反應速度。這些反應處理功能可在批反應模擬或一維對流、彌散、反應型運移模擬中使用。

2.GEM

GEM v.2009.13（Nghiem et al.，2004）是一款用來模擬利用CO₂和酸性氣體提高石油採收率的模擬器，該模擬器完全耦合了地球化學組成狀態方程。GEM採用一步求解法進行狀態方程的求解。GEM可以用來模擬：對流和彌散流體、油（或超臨界CO₂）、氣和鹹水間的平衡、水相物種間的化學平衡，以及礦物的動態溶解和沉澱。該模擬器採用自適應的隱式離散技術利用一維、二維或者三維模型來模擬孔隙介質中溶質的運移。油相和氣相用一個狀態方程來模擬，氣體在水相的溶解度採用亨利定律模型來計算。水通過蒸發進入到氣相、蓋層的穿透、熱效應和裂隙的封閉作用也可以利用GEM來模擬。

3.ECLIPSE

ECLIPSE是一個並行化的可以模擬黑油、組分、熱采等問題的成熟軟體。1994年，勝利石油管理局引進了ECLIPSE油藏數值模擬串列軟體，廣泛開展了從油藏到氣藏，從常用油田到特殊油氣田、從常規模擬研究到特殊模擬研究等多方面的應用。主要模塊有主模型、黑油、組分、熱采、流線法、運行平台和ECLIPSE Office等。

ECLIPSE是一個商業軟體，在使用中其內核部分是封閉的，使用者只能將其作為一個「黑箱」來操作。其不足之處有：不可能免費的獲得和隨意地使用和修改；無法耦合最前沿的地質流體熱力學模型；無法加入更多影響因素來研究具體問題。因此，ECLIPSE不適宜用於前沿科學研究。

4.TOUGH2/TOUGHREACT

TOUGH2是Transport of Unsaturated Groundwater and Heat（非飽和地下水流及熱流傳輸）的英文縮寫，是一個模擬一維、二維和三維孔隙或裂隙介質中，多相流、多組分及非等溫的水流及熱量運移的數值模擬程序。TOUGH2使用積分有限差分（Integral Finite Differences，IFD）（圖3-8）的方法來解決多相流動和多組分化學運移模擬中的空間離散化問題（Pruess et al.，1999s；Xu et al.，2004）。為了滿足大規模計算需要，Zhang et al.（2008）開發了TOUGH2的平行計算版本，即TOUGH2-MP。

該方法在對地質介質的離散化上是比較靈活的，允許使用不規則的網格，十分適合對多區域非均質系統和裂隙岩石系統中流體流動、運移和水岩相互作用的模擬。而對於規則的網格剖分，積分有限差分方法相當於傳統的有限差分法。其中，對於任意區域Vn，它的質量（對於水、氣體和其他化學組分）和能量（對於熱）守恆方程可以用積分的方式（式3-5）表達：

圖3-8 積分有限差分法中的空間離散化和幾何參數數據構成圖

中國二氧化碳地質封存選址指南研究

式中：下角標n為表示一個單元格；下角標m為表示和單元格n互相連接的網格m；Δt為時間步長；M_n為單元格n的平均質量或能量密度；A_nm為單元網格n和m交界的面段；F_nm為通過面積為A_nm的質量或能量通量；q_n為單元格n內單位體積的平均源匯率。

許天福等（1998）在TOUGH2的框架基礎之上，增加了多組分溶質運移和地球化學反應的模擬功能，形成了一套較為完善的可變飽和地質介質中非等溫多相流體反應地球化學運移模擬軟體——TOUGHREACT。該軟體不僅包括了TOUGH2的全部功能，而且適用於不同溫度、壓力、水飽和度、離子強度、pH值和氧化還原電位（Eh）等水文地質和地球化學條件下的熱-物理-化學過程。還可以應用於一維、二維或三維非均質（物理和化學的）孔隙或裂隙介質中的相關數值模擬研究。在理論上可以容納任意數量的以固相、液相或氣相存在的化學組分（但是在實際模擬中會受到計算能力和計算時間等硬體條件的限制），並且考慮了一系列化學平衡反應，如溶液中的配合反應、氣體的溶解或脫溶、離子吸附作用、陽離子交換及受平衡控制或反應動力學控制的礦物溶解或沉澱反應等。可以說TOUGHREACT、是TOUGH2的升級版，近年來在世界范圍內CO₂地質封存研究和工程實踐中得到了廣泛的應用。

除包含TOUGH2所有的功能外，TOUGHREACT還可以應用於一系列的反應性流體和地球化學遷移問題。比如：(1)伴隨Kd線性吸附和放射性衰變的污染物遷移問題；(2)在周圍環境條件下，自然界中地下水的化學演變；(3)核廢料處置地點評估；(4)深部岩層的沉積成岩作用；(5)CO₂地質處置。多相流體運動，多組分反應地球化學，各種封存形式封存量以及隨時間空間變化；(6)礦物沉積（如表生銅礦富集）；(7)自然和補給環境下熱水系統中的礦物變化。

通過最近幾年相關研究者的不懈努力，TOUGHREACT在實際應用中得到了進一步完善和提高，增加了部分新功能，如水相內部反應動力學和生物降解作用，改進了礦物-水反應表面積計算方法，以及氣-水反應中氣的活度系數的修正等。

5.PetroMod

由德國IES（Integrated Exploration System）公司研究開發的PetroMod多組分、多相態的多維含油氣系統模擬軟體綜合平台已被世界石油業所公認。該軟體融入了斷層活動性、鹽丘上涌和刺穿、火山岩的侵入、氣體擴散效應、油氣水三相運移和油氣吸附模型等相關技術。

該模擬軟體平台推出和採用的油氣運移組合模擬演算法（Hybird）是當今最先進的油氣運移模擬演算法，既可以保證模擬的精度，又可以極大地提高模擬的運算速度。其中的PetroFlow3D用於油氣運移、聚集、圈閉和散失等情況的模擬，同時PetroCharge Express為我們提供了基於圖件的油氣運移和圈閉模擬的快速分析工具。

6.MUFTE-UG

MUFTE-UG是MUFTE和UG.MUFTE的結合。MUFTE即多相流（Muliphase Flow）、運移（Transport）和能量（Energy）模型。該軟體包主要包括物理模型概念和孔隙裂隙介質中等溫和非等溫多相多組分流動和運移過程的離散方法（Helmig，1997；Helmig et al.，1998）。它能對裂隙孔隙介質進行離散性描述（Dietrich et al.，2005）。UG是非結構性網格（Unstructured Grid）的縮寫，它提供的數據結構能快速解算以平行、自適應多網格法為基礎的離散型偏微分方程。具有模塊化結構的MUFTE-UG很容易解決各種有特殊要求的問題。

模塊化結構的MUFTE-UG具有許多不同的環境與技術應用。例如，在環境應用領域，MUFTE-UG能夠模擬如下兩個問題。

（1）NAPL（非液相流體）向飽和與非飽和土壤的滲流。優化改進的修復技術在MUFTE中具有廣泛的研究和發展空間。

（2）地下CO₂的消散。CO₂以高溫高壓灌注地表以下幾百米的地層中，MUFTE-UG可用於非均質含水層（對流和彌散運移）中羽狀體演化評價，伴隨溫度效應（由於膨脹和壓縮）和組分間相互溶解（鹵水和CO₂）。

7.NUFT

NUFT（Nonisothermal Unsaturated-Saturated Flowand Transport model）是一套用來解決在多孔介質中多相、多組分非等溫流動和溶質運移過程中地下污染物運移的數值解法器。此軟體利用簡單的代碼來利用通用的實用程序和輸入文件的格式。最近，此代碼在Unix和DOS系統下運行成功。

該程序利用一套完整的有限差分空間離散法求解平衡方程組。每一個時間步長內利用Newton-Raphson方法求解非線性方程組，而在每一步迭代過程中利用直接解法和預共軛梯度法求解線性方程組。該模型可以解決一、二和三維水流及溶質運移問題。將來該模型會耦合進毛細滯後、非正交網格離散、有限單元剖分和固體非線性等溫吸附等功能。

（三）研究方法

通常情況下，CO₂地質封存數值模擬包括以下主要過程。

（1）建立概念模型：根據各種方法獲取的實際資料來概化和建立CO₂地質封存概念模型，包括邊界范圍、地層或儲蓋層高程、儲蓋層確定、參數及分區、源匯項、主要物理化學過程以及模型維度（一維、二維和三維）。

（2）建立數學模型：建立一套描述深部鹹水層中多相流動和多組分反應性溶質運移的偏微分方程組，包括初始條件和邊界條件問題。

（3）模型離散化：把概念模型中的各種信息通過網格剖分進行離散，形成大量的網格單元，然後通過有限差分、有限單元和積分有限差分等方法轉化成單元的質量和能量守恆方程組，再用多種方法將非線性方程組線性化，形成線性代數方程組，然後求解方程組。

（4）模型識別和校正：根據模型計算結果和實際監測數據進行對比擬合，適度合理調整參數，使模型能夠綜合反映實際情況。在歷史擬合過程中出現較大誤差，應重新檢查概念模型，修正概念模型。對所建模型進行參數敏感性分析，對於較敏感的參數應該慎重選取，甚至需要做大量的試驗來確定。

（5）模型預測：建立了可靠的模型後，便可以進行模擬預測。

數值模擬的關鍵是地質模型概化、計算精度和計算速度。由於計算的精度取決於離散的程度，而離散的程度又決定了計算的速度，這是一對矛盾，要根據解決問題的需要來選擇離散化的程度和計算速度。

CO₂在儲層中的運移、溶解以及與圍岩的化學反應形成了一個多相、多組分的反應體系，涉及的主要數學方程有超臨界CO₂-水的兩相流體運動控制方程、溶質運移控制方程和化學反應方程等。建立數值模型時，通常採用有限差分法、有限元法和積分有限差分法等。

由於實際應用時多採用已有的數值模擬軟體對CO₂地質封存的某一過程進行模擬，不涉及軟體的開發及程序代碼的編寫，只需根據研究的需要選擇合適的軟體進行模擬預測，而軟體一旦選定，數學模型和數值模型基本上已經確定。

『肆』 數值模擬的計算機方法

有限差分方法(FDM)是計算機數值模擬最早採用的方法，至今仍被廣泛運用。該方法將 求解域劃分為差分網格，用有限個網格節點代替連續的求解域。有限差分法以Taylor級 數展開等方法，把控制方程中的導數用網格節點上的函數值的差商代替進行離散，從而 建立以網格節點上的值為未知數的代數方程組。該方法是一種直接將微分問題變為代數 問題的近似數值解法，數學概念直觀，表達簡單，是發展較早且比較成熟的數值方法。 對於有限差分格式，從格式的精度來劃分，有一階格式、二階格式和高階格式。從差分 的空間形式來考慮，可分為中心格式和逆風格式。考慮時間因子的影響，差分格式還可 以分為顯格式、隱格式、顯隱交替格式等。目前常見的差分格式，主要是上述幾種形式 的組合，不同的組合構成不同的差分格式。差分方法主要適用於有結構網格，網格的步 長一般根據實際地形的情況和柯朗穩定條件來決定。

構造差分的方法有多種形式，目前主要採用的是泰勒級數展開方法。其基本的差分表達 式主要有三種形式：一階向前差分、一階向後差分、一階中心差分和二階中心差分等， 其中前兩種格式為一階計算精度，後兩種格式為二階計算精度。通過對時間和空間這幾 種不同差分格式的組合，可以組合成不同的差分計算格式。

有限元方法的基礎是變分原理和加權餘量法，其基本求解思想是把計算域劃分為有限個互不重疊的單元，在每個單元內，選擇一些合適的節點作為求解函數的插值點，將微分 方程中的變數改寫成由各變數或其導數的節點值與所選用的插值函數組成的線性表達式 ，藉助於變分原理或加權餘量法，將微分方程離散求解。採用不同的權函數和插值函數形式，便構成不同的有限元方法。有限元方法最早應用於結構力學，後來隨著計算機的發展慢慢用於流體力學的數值模擬。在有限元方法中，把計算域離散剖分為有限個互不重疊且相互連接的單元，在每個單元內選擇基函數，用單元基函數的線形組合來逼近單元中的真解，整個計算域上總體的基函數可以看為由每個單元基函數組成的，則整個計算域內的解可以看作是由所有單元上的近似解構成。在河道數值模擬中，常見的有限元計算方法是由變分法和加權餘量法發展而來的里茲法和伽遼金法、最小二乘法等。根據所採用的權函數和插值函數的不同，有限元方法也分為多種計算格式。從權函數的選擇來說，有配置法、矩量法、最小二乘法和伽遼金法，從計算單元網格的形狀來劃分，有三角形網格、四邊形網格和多邊形 網格，從插值函數的精度來劃分，又分為線性插值函數和高次插值函數等。不同的組合 同樣構成不同的有限元計算格式。對於權函數，伽遼金(Galerkin)法是將權函數取為逼近函數中的基函數 ；最小二乘法是令權函數等於餘量本身，而內積的極小值則為對代求系數的平方誤差最小；在配置法中，先在計算域 內選取N個配置點 。令近似解在選定的N個配置點上嚴格滿足微分方程，即在配置點上令方程餘量為0。插值函數一般由不同次冪的多項式組成，但也有採用三角函數或指數函數組成的乘積表示，但最常用的多項式插值函數。有限元插值函數分為兩大類，一類只要求插值多項式本身在插值點取已知值，稱為拉格朗日(Lagrange)多項式插值；另一種不僅要求插值多項式本身，還要求它的導數值在插值點取已知值，稱為哈密特(Hermite)多項式插值。單元坐標有笛卡爾直角坐標系和無因次自然坐標，有對稱和不對稱等。常採用的無因次坐標是一種局部坐標系，它的定義取決於單元的幾何形狀，一維看作長度比，二維看作面積比，三維看作體積比。在二維有限元中，三角形單元應用的最早，近來四邊形等參元的應用也越來越廣。對於二維三角形和四邊形電源單元，常採用的插值函數為有Lagrange插值直角坐標系中的線性插值函數及二階或更高階插值函數、面積坐標系中的線性插值函數、二階或更高階插值函數等。

對於有限元方法，其基本思路和解題步驟可歸納為
(1)建立積分方程，根據變分原理或方程餘量與權函數正交化原理，建立與微分方程初邊值問題等價的積分表達式，這是有限元法的出發點。
(2)區域單元剖分，根據求解區域的形狀及實際問題的物理特點，將區域剖分為若干相互連接、不重疊的單元。區域單元劃分是採用有限元方法的前期准備工作，這部分工作量比較大，除了給計算單元和節點進行編號和確定相互之間的關系之外，還要表示節點的位置坐標，同時還需要列出自然邊界和本質邊界的節點序號和相應的邊界值。
(3)確定單元基函數，根據單元中節點數目及對近似解精度的要求，選擇滿足一定插值條 件的插值函數作為單元基函數。有限元方法中的基函數是在單元中選取的，由於各單元 具有規則的幾何形狀，在選取基函數時可遵循一定的法則。
(4)單元分析：將各個單元中的求解函數用單元基函數的線性組合表達式進行逼近；再將 近似函數代入積分方程，並對單元區域進行積分，可獲得含有待定系數(即單元中各節點 的參數值)的代數方程組，稱為單元有限元方程。
(5)總體合成：在得出單元有限元方程之後，將區域中所有單元有限元方程按一定法則進 行累加，形成總體有限元方程。
(6)邊界條件的處理：一般邊界條件有三種形式，分為本質邊界條件(狄里克雷邊界條件 )、自然邊界條件(黎曼邊界條件)、混合邊界條件(柯西邊界條件)。對於自然邊界條件， 一般在積分表達式中可自動得到滿足。對於本質邊界條件和混合邊界條件，需按一定法 則對總體有限元方程進行修正滿足。
(7)解有限元方程：根據邊界條件修正的總體有限元方程組，是含所有待定未知量的封閉 方程組，採用適當的數值計算方法求解，可求得各節點的函數值。

有限體積法（Finite Volume Method）又稱為控制體積法。其基本思路是：將計算區域劃分為一系列不重復的控制體積，並使每個網格點周圍有一個控制體積；將待解的微分方程對每一個控制體積積分，便得出一組離散方程。其中的未知數是網格點上的因變數的數值。為了求出控制體積的積分，必須假定值在網格點之間的變化規律，即假設值的分段的分布的分布剖面。從積分區域的選取方法看來，有限體積法屬於加權剩餘法中的子區域法；從未知解的近似方法看來，有限體積法屬於採用局部近似的離散方法。簡言之，子區域法屬於有限體積發的基本方法。

有限體積法的基本思路易於理解，並能得出直接的物理解釋。離散方程的物理意義，就 是因變數在有限大小的控制體積中的守恆原理，如同微分方程表示因變數在無限小的控 制體積中的守恆原理一樣。 限體積法得出的離散方程，要求因變數的積分守恆對任意一組控制體積都得到滿足，對整個計算區域，自然也得到滿足。這是有限體積法吸引人的優點。有一些離散方法，例如有限差分法，僅當網格極其細密時，離散方程才滿足積分守恆；而有限體積法即使在粗網格情況下，也顯示出准確的積分守恆。就離散方法而言，有限體積法可視作有限單元法和有限差分法的中間物。有限單元法必須假定值在網格點之間的變化規律（既插值函數），並將其作為近似解。有限差分法只考慮網格點上的數值而不考慮值在網格點之間如何變化。有限體積法只尋求的結點值，這與有限差分法相類似；但有限體積法在尋求控制體積的積分時，必須假定值在網格點之間的分布，這又與有限單元法相類似。在有限體積法中，插值函數只用於計算控制 體積的積分，得出離散方程之後，便可忘掉插值函數；如果需要的話，可以對微分方程 中不同的項採取不同的插值函數。

『伍』 數值模擬流程

不同的軟體進行數值模擬時所需的參數、計算方法、剖分格式等不盡相同，數值模擬的過程也不同，但大致相同，本文以TOUGHREACT為例介紹CO₂地質儲存數值模擬的流程。

（一）研究范圍的確定

一般情況下，獨立的天然水文地質系統是計算區最好的選擇，它具有自然邊界，便於較准確地利用其真實的邊界條件，避免人為邊界在資料提供上的困難和誤差。但是在實際工作中，常常不能完全利用自然邊界，這就需要充分利用勘察和長期觀測資料等建立人為邊界。在確定計算區域時，除了保證范圍足夠大以外，還應使假定的邊界條件盡可能接近真實狀態。

計算范圍的劃定應充分考慮研究目的、區域地質構造、儲層岩性、儲層岩石礦物組成及地下水化學成分等多方面因素。數值模擬時間根據研究目的不同稿謹具有不同的時間尺度。就CO₂地質儲存數值模擬而言，如果不考慮地球化學作用，封存系統在1000年數量級的模擬時間內基本上已達到平衡或穩定。在劃定邊界時還應考慮CO₂在儲層中的擴散距離，與研究區地質模型的孔隙度，滲透率等參數關系密切。為了保證所選模型範圍邊界在模擬期內不影響模擬結果，盡量通過具有相同地質條件的天然CO₂氣田（藏）進行類比，確定大體的計算范圍的邊界。如果考慮地球化學反應，由於CO₂注入引發的水－岩－氣反應對圍岩岩性改變較顯著，制約著CO₂注人的速度和徑向運移的距離等。

（二）明確研究目的

在進行數值模擬以前首先要明確利用數值模擬技術要解決什麼樣的問題。對於CO₂地質儲存工程而言，進行數值模擬的目的主要是在CO₂地質儲存工程實施前，通過數值模擬技術對工程的選址、方案設計進行優化，工程實施期技術指導、運沖森行期監測及後期CO₂泄漏的風險評估等進行預測，以指導項目科學、合理地實施，將CO₂泄漏風險降至最低。

研究目的決定著前期資料的收集類型、地質建模的側重點、地質模型離散的精密度以及初始、邊界條件的處理方式等過程。

（三）資料的收集整理

1）通過遙感、綜合地質調查、物探、鑽探和各類樣品測試分析等手段獲取場地深部地層岩性、地質構造、水文地質、水文地球化學、岩石礦物資料和數據；

2）搜集和分析CO₂地質儲存場地地質岩性、區域構造格架、活動斷層與地震活動情況等；

3）採用鑽井岩心、測井和地震反射方法，調查CO₂地質儲存場地目標儲層和蓋層的空間分布形態，埋深、厚度和規模等；

4）使用X射線衍射、掃描電鏡等方法研究分析封存場地岩石礦物組成、孔隙結構特徵及其物理化學性質；

5）通過採取淺部、深部含水層水樣進行水質全分析，獲得儲蓋層地層水及淺部含水層初始水化學成分。

不同的數值模擬軟體其數學模型的數值解法不同，空間離散方式也不盡相同，所需的模型參數也有一定的差異，表9-1即為TOUGHREACT數值模擬所需要的主要參數。

表9-1 CO₂地質儲存模擬過程中需要的主要參數（以TOUGHREACT為例）

圖9-3 網格剖分

網格剖分對計算的精度及計算的效率有很重要的影響。精度越高對模擬結果刻畫的越精細，但是數據的計算量越大，對計算機的要求也越高。建議在進行地質模型剖分時先採用較粗的網格剖分，如果這種剖分方式下模擬結果合理然後再進行精細化剖分，用於對模擬結果更加詳細的刻畫。

2.參數和初始條件

初始條件是指在初始時刻（t=0）時研究區內求解數學模型主要狀態變數的初始值。選擇的應用軟體不同所需的狀態變數數量、種類不同。如TOUGHREACT所需的初始主要狀態變數包括壓力、溫度和組分濃度的空間分布。地質參數包括孔隙度、滲透率、密度、壓力、溫度、毛細壓力等參數值。這些數值一部分採用室內實驗測得，另一部分採用參考文獻的經驗值；地層水的化學成分的初始值採用實際地層水的化學分析，主要是8大離子的濃度、鹽度和pH 等。如果研究區深部地層中的水樣難以獲得，如蓋層，則採用靜態平衡的方法，利用具有與儲層相同鹽度的鹹水與含有原生礦物的地層岩石在原地層環境下進行化學反應，獲取平衡狀態下的地層水化學成分的初始值；通過岩礦分析、電子掃描、Ⅹ衍射等手段，獲得組成CO₂地質儲層蓋層的原生礦物成分體積含量初始值，並根據原生礦物的組成合理判斷次生礦物。

從原則上講，初始時刻是可以任意取定的，只要該時刻所需的參數和狀態變數值已知即可。因此我們不應該把初始條件理解為研究系統的初始狀態。具體如何取，應該視問題的需要、資料來源、計算方便與否等因素而定。

3.邊界條件

邊界條件是某一實際問題數學模型具有定解的必要條件之一。地下水流問題和溶質運移問題邊界條件的定義不盡相同，但一般概化為以下三種。

（1）一類邊界條件（Dirichlet條件）

解決水流問題時，此類邊界條件為在邊界上所有點的水頭是給定的；對於溶質運移問題，一類邊界條件是指研究區邊界上的溶質濃度分布已知。解決CO₂—水兩相流動問題時，此類邊界條件為在邊界上所有點的壓力是給定的。

（2）二類邊界條件（Neumann條件）

當已知某一邊界的單位面積流入或流出的流量時，可視作解決流動問題的二類邊界；相對溶質運移來講，此類邊界又稱給定彌散通量邊界，即邊界上的彌散通量隨時間變化規律已知。

（3）三類邊界（Cauchy條件）

當研究區一部分滿足一類Dirichlet條件，而另一部分滿足二類Neumann條件時，這類問題稱為混合邊界問題，稱為三類邊界。對溶質運移而言，此類邊界為邊界上溶質通量隨時間變化規律已知。

在CO₂地質儲存數值模擬過程中，由於儲層地層多在800m以下，地質模型的頂部和底部根據實際需要可以處理為不透邊界；為了避免邊界對模擬結果的影響，研究區的范圍一般比實際CO₂所能運移到的范圍大得多，因此，在處理四周邊界時一般設置為無窮一類邊界或不透邊界。在確定邊界條件時，應根據水文地質條件以及現有的資料來綜合考慮。

4.源匯項處理

在多孔介質中流動和溶質運移的問題中，對流、水動力彌散和溶質源或／和匯，是決定含水層中任一內點上溶質質量時變率的兩大因素。源匯項問題在水質與水量計算中以及正確處理對流－彌散方程和滲流基本微分方程中佔有重要地位。作為源匯項的方式很多，如越流補給、含水層彈性釋放補給以及抽（注）井的補給等。

對於深部鹹水層CO₂地質儲存系統而言，系統的頂部一般為具有低滲、低孔的泥岩、頁岩等緻密性岩層，越流補給較難發生。整個CO₂地質儲庫系統的源匯項主要指對流（如側向邊界）和抽（注）井。

（八）模型的校正與驗證

模型識別是建立地下流體數值模型最重要的環節之一，正確理解和進行擬合對於提高數值模型的模擬性是至關重要的。在有實測結果的情況下如示範工程，可將模擬結果與實測結果進行比較，對相關參數進行適當合理的調整，使模擬結果在給定的誤差范圍內與實測結果吻合。若誤差較大，應該重新檢驗概念模型的可靠性，甚至重新建立概念模型。在識別校正以後，應採用校正好的模型繼續計算，並與未用來識別校正的實際數據比較，驗證模型的准確性和可靠性。若存在較大誤差，需重復前面的過程。在沒有實測結果的情況下，數值模型的可靠性可通過類比相關資料或根據個人經驗和理論判斷。

（九）模擬預測

模型預測是實施數值模擬技術的主要目的。對於CO₂地質儲存工程而言，由於CO₂地質儲存技術的提出為時尚短，針對CO₂在深部鹹水層中的運移、擴散、與地層水和圍岩產生的化學反應，以及由於CO₂灌注引起的儲蓋層物理、化學性質變化研究均處於研究和發展階段。因此，在工程實施過程中急需具有技術指導性的工具產生，避免造成投資的浪費和CO₂泄漏等風險的出現。

利用經過識別校正與驗證過的數值模型對CO₂地質儲存過程進行模擬預測，有針對性地對模擬數據進行後期處理，如統計分析、比較等手段對結果進行解譯，以此達到場地的優選，目標儲層灌注能力、儲存潛力的評估，CO：擴散運移途徑和速度、不同捕集方式封存量及它們之間的時空轉化等過程的詳細刻畫與模擬模擬等目的。同時可以預測CO₂在已有、重新激活或新生成的裂隙中逃逸的可能性及時間、CO₂泄漏風險評估以及評價CO₂泄漏對淺層地下水的水質、水量及對地表環境的影響等。

上述結果的分析只是數值模擬技術所能解決問題的冰山一角。對於數值模擬結果的處理要根據所研究的目的進行有針對性的提取和解譯。通過對處理後的數據進行總結分析，發現問題從而解決問題，並掌握內在規律，為CO₂地質儲存工程的前期設計、工程實施、中期監測管理提供理論支持和科學的技術指導，並可以提前開展風險預測，盡早制定預案防範CO₂地質儲存工程實施及運行過程中可能出現的隱患。

『陸』 模擬方法簡介

4.4.1.1 FLAC^3D的基本原理

數值計算一直是工程地質學分析方法中的一個重要手段，特別是近年來，隨著計算機技術和數值計算方法的迅速發展，從二維到三維，從靜態到動態，計算精度、可靠度不斷提高，加之大量工程實例的驗證，使數值模擬計算方法已日益成為岩體穩定評價中不可或缺的重要方法之一。本研究利用國際工程地質界公認的最新通用軟體FLAC^3D2.0作為工具，對典型岩溶塌陷類型的工程環境效應進行評價。

FLAC^3D分析在求解中使用如下3種計算方法：

（1）離散模型方法，連續介質被離散為若干六面體單元，作用力均被集中在節點上。

（2）有限差分方法，變數關於空間和時間的一階導數均用有限差分來近似表達。

（3）動態鬆弛方法，由質點運動方程求解，通過阻尼使系統運動方程衰減至平衡狀態。

4.4.1.1.1 空間導數的有限差分

快速拉格朗日分析採用混合離散方法，將區域離散為常應變六面體單元的集合體，又將每個六面體看作以六面體角點為角點的常應變四面體的集合體，應力、應變、節點不平衡力等變數均在四面體上進行計算，六面體單元的應力應變取值為其內四面體的體積加權平均，這種方法既避免了常應變六面體單元常會遇到的位移剪切鎖死現象，又使得四面體單元的位移模式可以充分適應一些本構的要求，如一四面體，節點編號為1到4，第n面表示與節點n相對的面，設其內一點的速率分量為υ_i，由高斯公式得

岩溶塌陷機理及其預測與評價研究

式中：V為四面體的體積；S為四面體的外表面；n_j為外表面的單位法向向量分量。對於常應變單元，υ_i為線性分布，n_j在每個面上為常量，由式（4-1）可得：

岩溶塌陷機理及其預測與評價研究

式中：上標l表示節點l的變數；上標（l）表示面l的變數。

4.4.1.1.2 運動方程

快速拉格朗日分析以節點為計算對象，在時域內求解，節點運動方程可表示如下：

岩溶塌陷機理及其預測與評價研究

式中

（t）為在t時刻l節點在i方向的不平衡力分量，可由虛功原理導出；m^l為l節點的集中質量。對於靜態問題，採用虛擬質量以保證數值穩定，而對於動態問題則採用實際的集中質量。

將上式左端用中心差分來近似，則可得：

岩溶塌陷機理及其預測與評價研究

4.4.1.1.3 應變、應力及節點不平衡力

快速拉格朗日分析由速率來求某一時步的單元應變增量，即

岩溶塌陷機理及其預測與評價研究

有了應變增量，即可由本構方程求出應力增量，進而得到總應力。

4.4.1.1.4 阻尼力

對於靜態問題，在式（4-3）的不平衡力中加入了非粘性阻尼，以使系統的振動逐漸衰減直至達到平衡狀態（即不平衡力接近零），此時式（4-3）變為：

岩溶塌陷機理及其預測與評價研究

阻尼為：

岩溶塌陷機理及其預測與評價研究

式中：α為阻尼系數，其默認值為0.8；而符號函數

岩溶塌陷機理及其預測與評價研究

4.4.1.1.5 計算循環

由以上可以看出快速拉格朗日分析的計算循環如圖4-2所示。

圖4-2 FLAC^3D軟體計算循環圖

4.4.1.2 FLAC^3D的功能和應用范圍

盡管FLAC^3D的計算公式源於有限差分方法，但其計算結果與有限元方法的計算結果（對於常應變四面體）相同，而且它與現行的數值方法相比有著明顯的優點（黃潤秋等，1991、1995、1999a、1999b），其特點包括：

（1）「混合離散化」（mixed discretization）技術的應用，更能精確和有效地模擬計算材料的塑性破壞和塑性流動，在力學上比常規有限元的數值積分更為合理。

（2）全部使用動力運動方程，即使在模擬靜態問題時也如此。因此，它可以較好地模擬系統的力學不平衡到平衡的全過程，實現動態的模擬過程。

（3）求解中採用「顯式」差分方法大大節約了計算時間，特別對求解任意的非線性應力-應變問題尤為重要。同時它不需要存儲較大的剛度矩陣，因此，它與一般的差分分析方法相比，既節約了計算機的內存空間又減少了運算時間，大大提高了解決問題的速度。

（4）物體由多面體單元所表示，可以通過調整三維網格的方法，以適應研究體真實的形狀，每個單元力學行為是對應力-應變法則和邊界力、約束條件的響應。材料能產生屈服和流動，而且網格也能變形（in large-strain mode），並隨材料移動。

（5）強大的後處理功能。能根據需要輸出設定階段的應力應變成果，提取各工程部位的應力、應變的量值，同時還可以根據追蹤功能（history）獲取所需變數值的歷時變化曲線。

研究中使用的FLAC^3D所使用的模型主要包括近10種模型，常用的有以下幾種：

（1）空模型（null）。用空材料來表示開挖後的材料，在空單元中的應力被自動地賦零。

即：

岩溶塌陷機理及其預測與評價研究

用此種模型模擬實際中的洞穴。

（2）彈性均質模型（elastic，isotropic）。在FLAC^3D中，對於基岩主要採取了這種模型，而對於非均質的介質再根據其不均勻情況進行特殊的定義。在這種模型中，據虎克（Hooke）定律，有：

Δσ_ij=2GΔε_ij+α₂Δε_kkδ_ij

在此Einstein求和中，δ_ij為位移，α₂為與材料的體積模量及剪切模量相關的材料系數：

岩溶塌陷機理及其預測與評價研究

新的應力值通過下式獲得：

岩溶塌陷機理及其預測與評價研究

（3）彈性非均質模型（elastic，orthotropic），主要用於非均質情況的模擬。

（4）Drucker-Prager塑性模型。在FLAC^3D中此模型使用的破壞准則主要為Drucker-Prager拉張破壞准則，用以評價計算中的塑性破壞現象。

（5）摩爾-庫侖模型（Mohr-Coulomnb）。在整個數值模擬過程中，對於土層、砂土、粉土等均採用了摩爾-庫侖模型作為計算模型，即將其視為彈塑性材料來解決。在此模型中，以Mohr-Coulomnb准則作為材料的破壞准則。